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ABSTRACT 



In this thesis different numerical methods, as well as applications of the 
methods to a number of current problems in relativistic astrophysics, are 
presented. 

In the first part the theoretical foundation and numerical implementation 
of a new general relativistic magnetohydrodynamics code is discussed. A new 
form of the equations of motion using global coordinates, but evolving the 
dynamical variables from the point of view of a local observer is presented. 
No assumptions are made about the background metric and the design is 
ready to be coupled with methods solving the full Einstein equations. 

In the second part of the thesis important results concerning the under- 
standing of coUisionless shocks, obtained from experiments with a relativistic 
charged particle code, are presented. Relativistic coUisionless shocks arc im- 
portant in a range of astrophysical objects; in particular in gamma ray burst 
afterglows and other relativistic jets. It is shown that a strong small scale, 
fluctuating, and predominantly transversal magnetic field is unavoidably gen- 
erated by a two-stream instability. The magnetic energy density reaches a 
few percent of equipartition. 

A new acceleration mechanism for electrons in ion-electron coUisionless 
shocks is proposed. The mechanism is capable of creating a powerlaw electron 
distribution in a coUisionless shocked region. The non-thermal acceleration 
of the electrons is directly related to the ion current channels generated by 
the two-stream instability and is local in nature. Thus the observed radiation 
field may be tied directly to the local conditions of the plasma and could be 
a strong handle on the physical processes. 

Experiments of colliding pair plasmas are presented and the formation of 
a macrophysical shock structure is observed. A comparable relativistic fluid 
simulation is performed and good agreement is found, implying that the full 
structure of the shock has been resolved. The extent of the shock transition 
region in a pair plasma is estimated to 50-100 electron skin depths. 

In the third part of the thesis a new particle-in-cell code is discussed. It 
solves the full Maxwell equations, together with direct microphysical particle- 
particle interactions, such as relativistic scattering, pair production, decay, 
and annihilation of particles. The inclusion of such relativistic interaction 
processes makes it possible to extract self consistent synthetic photon spec- 
tra directly from the numerical experiments, thereby gaining the ability to 
directly compare models with observations. 
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1. INTRODUCTION & OVERVIEW 



During the last decade we have seen fundamental advances in the observation 
of compact objects, active galactic nuclei, gamma ray bursts, and other ob- 
jects characterised by their extreme physical conditions and emittence of light 
over the full electromagnetic spectrum. This branch of astrophysics has aptly 
been named "extreme astrophysics" , and advances in the field are driven by 
the technical development and launch of new satellites, such as Bcppo/Sax, 
Chandra, XMM and Swift, and the construction of powerful ground based 
facilities, such as the HESS telescope and the Auger observatory to measure 
X-rays and gamma rays. Moreover the technique of combining radio tele- 
scopes to perform interferometric observations with synthetic dishes compa- 
rable to the entire globe has played an important role for resolving the inner 
engines of active galactic nuclei (AGN). 

In 1997 the first afterglow from a gamma ray burst (GRB) was observed, 
placing GRBs firmly out of reach at cosmological distances and earning them 
the title as the most violent explosions in the Universe. Very high energy 
gamma rays have also been observed from AGNs, and with the increasing 
resolution of high frequency radio interferometers, we will soon be able to 
resolve the launching region of the jets associated with AGNs, only a few 
Schwarzchild radii from the central supermassive black hole. 

In the decade to come we can foresee that two entirely new windows to the 
Universe will be opened to complement the observations of electromagnetic 
radiation and cosmic rays that arc made today: On the south pole the Ice- 
Cube project will detect cosmic neutrinos generated in the core of supernovae 
and possibly in GRBs and other cataclysmic events, while laser interferom- 
eters on the ground, such as LIGO, VIRGO and GEO 600, together with 
the space interferometer LISA, will have reached levels of sensitivity where 
the gravitational waves from the coalescence of compact objects and super 
massive black holes may be detected. 

A decade ago cosmology was still the branch of astrophysics where one 
could get along with back-of-the-envelope calculations, since fundamental 
parameters such as the Hubble expansion rate, the age and the matter con- 
tent of the Universe were all quoted with error bars as large as the numbers 
themselves. This is all in the past now. The Hubble space telescope has 
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finally determined the expansion rate. Observations of supernovae of type la 
at moderate and high redshifts have led to the surprising conclusion that the 
Universe is in fact accelerating in its expansion. The Boomerang and Max- 
ima balloon missions and later the WMAP telescope have nailed down fluc- 
tuations in the cosmic microwave background radiation (CMBR) with high 
precision and determined the overall geometry of the Universe to be. . . flat! 
Euclid was right. The pieces in the cosmological puzzle are slowly falling 
into place. Current and future dedicated facihties to observe the CMBR, 
together with large scale galaxy redshift surveys such as the SLOAN digital 
sky survey and the 2DF survey, will give strong limits on the distribution of 
matter and flelds in the early Universe. 

It is thus fair to say that both extreme astrophysics and cosmology, to- 
gether known as relativistic astrophysics, are in a golden age and are slowly 
but flrmly entering the realm of "messy astrophysics" , where predictions can- 
not be based on sketchy ideas anymore but instead detailed physical models 
must be worked out, tested, and validated or falsifled by observations. 

Parallel to the development in observational relativistic astrophysics, there 
has been a revolution in the tools employed by theoretical astrophysicists. 
The computational power has for decades been rising exponentially, dou- 
bling every 18 months in accordance with Moores law. At the end of the 
nineties three-dimensional computer models of astrophysical objects became 
affordable, and for some time computer modelling has been indispensable in 
understanding the Universe. 

In order to interpret observations, we have to develop theories that in 
simple terms grasp the central physical mechanisms and let us understand 
how fundamental parameters are related. As the observations become more 
complicated, and the quality of the data improves, so must the theories to be 
successful in explaining these new details. Astronomy is different from other 
natural sciences, in that we cannot perform experiments in the laboratory, 
and in most cases the timescales are so long that we cannot even wait and 
watch them unfold in the Universe. 

In compact objects and in the early Universe many different physical pro- 
cesses play important roles to shape the flnal picture, ranging from the large 
scale fluid dynamics, the curvature of space, the interactions in the plasma 
between matter and electromagnetic fields, all the way to the microphysi- 
cal generation and scattering of the light, which, ultimately, is observed on 
Earth. The computer gives us, as a complement to observations, the ability 
to create models, and in contrast to the real Universe, we can spin our models 
around and visualise the data in three dimensions, instead of the projected, 
two-dimensional view which is the only one that the real Universe offers. In 
this sense computer models have become the virtual laboratory of the astro- 
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physicist. The physical insights gained from these models are essential, and 
often the complexity of the phenomena leaves us at loss without access to 
such models. 

1.1 A Swiss Army Knife for Relativistic Astrophysics 

In this thesis I present the application, development and implementation of 
several computer codes which may be used to model relativistic astrophysics. 
They span a range of scales and interactions. The GrMHD code, presented 
in Chapter 2, may be used to describe the flow of matter from cosmological 
scales down to the scales of black holes. The charged particle code, used in 
Chapters 3-5, is applied to understanding the small scale structure in colli- 
sionless shocks. Finally, the photon plasma code, presented in Chapter 6, will 
enable us to study a fuller range of plasma physics, including microphysical 
interactions, scatterings and the detailed propagation of radiation. 

1.2 General Relativistic Magneto-Hydrodynamics 

In Chapter 2 I present a reformulation of the equations of motion for general 
relativistic magnetohydrodynamics (GrMHD) that is well suited for numeri- 
cal purposes, and the implementation in a three-dimensional numerical code 
that solves the equations. Before starting the implementation of the code, 
I carefully considered the approaches employed in the handful of existing 
codes worldwide. My main idea has been to make a conscious split between 
the reference frame in which we measure our coordinates, and the reference 
frame in which we measure our physical variables. 

The coordinate system, naturally, has to cover the whole physical domain. 
In the case of compact objects, it is normal to use a coordinate system 
connected to observers at infinity. But there is no a priori reason why we 
should measure physical variables, such as density, velocity, and internal 
energy, as seen by observers at infinity. If one measures them in a locally 
defined frame, which is related to a local inertial frame, then the physics, by 
the equivalence principle, becomes almost like the physics of special relativity, 
for arbitrary background space times. 

All equations have been derived without placing any constraints on the 
metric tensor. It is important to keep everything completely general, to allow 
the code in the future to be enhanced with procedures that solve the Einstein 
equations and evolve the metric tensor. 

The code is based on finite difference techniques, and to handle discon- 
tinuities we have to include some form of artificial viscosity to enhance the 
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Fig. 1.1: A relativistic jet. From left to right: The density, the pressure and the 
total four velocity. The jet head has partly destabilised and is creating a 
complex bubble of vortices in the cocoon of the jet. 



entropy across shock fronts. I have chosen to employ the correct physical de- 
scription of shear viscosity, to enforce energy and momentum conservation. 
The full shear viscosity is very complicated, and in the general case of an 
arbitrary space time it would be impractical if we did not use variables tied 
to a local frame of reference. 

The code has been subjected to an extensive testbed of hydrodynamic and 
magnetohydrodynamic problems, and it is demonstrated that it can handle 
large gradients in density, pressure and the magnetic field. As an example 
Fig. 1.2 shows a highly relativistic shock tube problem, with two shock waves 
travelling toward each other, involving a relative pressure jump of 10^. The 
solution is shown at different resolutions with the analytic solution overplot- 
ted. This test is described further as problem III in Chapter 2. Moreover, as 
an example of the capabilities the code, I use it in two relevant astrophysical 
applications, one of them being the injection of a thin and hot relativistic jet 
into an ambient medium shown in Fig. 1.1. 

I have implemented a fully three-dimensional version of the code. The 
code is parallelised, has been tested on several supercomputers, and has been 
shown to yield excellent performance on hundreds of CPUs. 

1.3 Magnetic Field Generation in CoUisionless Shocks 

Chapter 3 was published by Frederiksen, Hededal, Haugb0lle and Nordlund 
[26]. Using three-dimensional particle simulations we report on the evolution 
of an ion-electron dominated counter-streaming collisionless shock. Our ex- 
periment consists initially of two populations. An in-streaming population. 
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Fig. 1.2: A highly relativistic shock tube. The solution is shown at t = 0.2, just 
before the two shock waves collide. 
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Fig. 1.3: Left: The initial conditions for our experiment. Right: Electron (top) 
and ion (bottom) currents, averaged over the x-direction. The plasma is 
streaming from left to right. 



with a Lorentz boost of F = 3 upstream of the shock interface, and a pop- 
ulation at rest downstream of the shock interface (see Fig. 1.3 to the left). 
It is predicted theoretically, that colliding coUisionless plasmas are suscepti- 
ble to the Weibel- or two-stream instability. Microscopic fluctuations in the 
magnetic field deflect the particles that in turn enhance the fluctuation and 
an exponential growth of the fluctuations results in the generation of strong 
current channels. In our simulations this is confirmed and we observe the 
instability develop in the shock interface. In Fig. 1.3 to the right is shown 
the current densities at late times. Associated with the current channels 
is a strong transversal magnetic field. The magnetic field energy density 
reaches a few percent of the kinetic energy in the in-coming beam. For an 
ion-electron plasma this is in fact a two stage process. Initially when the 
electrons encounter the shock interface, being the lighter particles they are 
rapidly deflected into first caustic surfaces and then current channels. The 
magnetic field keeps growing in scale and strength, until the ions undergo 
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the same process and similarly ion current channels are formed. Because 
of charge separation, the electrons will be attracted to the ions, and the 
electron instabihty is quenched. Instead the electrons start to Debye shield 
the ions, forming a fuzzy cloud around the ion channels (see Fig. 1.3). The 
Debye shielding partly neutralises the ion channels, and helps stabilise the 
evolution. The electrons are fully thermalised, but the ions are only slightly 
deflected from their initial distribution, due to the strong shielding of the 
electrons, and thermalisation might be significantly slower than predicted 
simply by extrapolating with the mass ratio. The ion current channels mu- 
tually attract each other and a self similar merging process commences, where 
neighbouring channels merge to form larger channels. With the capacity of 
current computers, the ions cannot be followed all the way to thermalisation, 
and merging of current channels is ongoing when they reach the end of the 
box and stream out at the open boundary. 

To generate the radiation seen in observations of GRB afterglows, a mag- 
netic field containing 10~^ — 10~^ of the kinetic energy is required. The 
two-stream instability seen to occur in our experiments is a strong candidate 
for explaining the generation of this magnetic field, since it is unavoidable in 
coUisionless shocks with low degrees of magnetisation. It then follows that 
the magnetic field cannot be taken as a free parameter, but is a consequence 
of the parameters of the shock, such as the inflow velocity and the density 
contrast. These findings do not only pertain to GRB afterglows, but also 
imply that magnetic field generation may be an important ingredient in all 
weakly magnetised coUisionless shocks, and therefore occurs in a range of 
objects from supernovae remnants to internal shocks in outflows from AGN, 
all harbouring coUisionless shocks. 

1.4 Non-Fermi Powerlaw Acceleration in Astrophysical 

Plasma Shocks 

In Chapter 4 I present the results published by Hededal, Haugb0lle, Pred- 
eriksen and Nordlund [32]. We study highly relativistic charged ion-electron 
particle dynamics in coUisionless shocks. The numerical experiment reported 
on here is different from the one in Chapter 3 in that the in-streaming plasma 
has a higher Lorentz factor (F = 15) and the computational box employed 
is about 3 times longer in the streaming direction, enabling us to follow the 
process further downstream of the shock interface and for a longer period of 
time, until the shock structure has been more fully developed. 

We flnd a powerlaw distribution of accelerated electrons, which turns out 
to originate from an acceleration process that is a direct consequence of the 
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two-stream instability observed in Chapter 3 and is local in nature. The 
electrons are accelerated and decelerated when passing through the cores of 
the ion current channels generated by the two-stream instability, and the 
process is fundamentally different from recursive acceleration processes, such 
as Fermi acceleration. We find a powerlaw slope of 2.7, in concordance with 
that inferred from observations of the afterglow in gamma ray bursts, and 
the process may explain more generally the origin of part of the non-thermal 
radiation from relativistic jets, supernovae remnants and shocked inter- and 
circum-stellar regions. 

When two coUisionless plasmas interpenetrate, current channels are formed 
through the two-stream instability. The ion current channels dominate the 
dynamics, due to the heavier mass of the ions, and downstream of the shock 
the channels merge in a hierarchical manner forming increasingly stronger 
patterns. The electrons act to Debye shield the channels yielding charge 
neutrality at large distances. At distances less than the Debye length the ion 
channels are surrounded by an intense transverse electric field that accelerate 
the electrons toward the channels and then decelerate them, when they move 
away from the channel. This can be seen in Fig. 1.4, where we in part (A) 




Fig. 1.4; (A) Ray traced electron paths (red) and current density (blue). The 
colours of the electron paths reflect their four velocity according to the 
colour table in the inset (B). The shadows are equivalent to the x and 
y projections of their paths. The ion current density is shown with blue 
colours according to the colour table in the inset. The inset also shows 
the ion current density (blue) integrated along the x axis with the spatial 
distribution of fast moving electrons (red) over plotted. 

have ray traced two selected electron paths and colour coded them according 
to the velocity and in part (B) have shown the spatial distribution on the 



1.4 Non-Fermi Powerlaw Acceleration in Astrophysical Plasma Shocks 



8 



fastest electrons in the box overplotted on top of the ion current distribu- 
tion. Notice the strong correlation between fast moving electrons and high 
ion current density. 

To analyse the process quantitatively we have constructed a toy model, 
idealising the ion channel as a solid cylinder of moving ions. Given an elec- 
tron we can calculate the maximal energy gained in the acceleration towards 
the cylinder (see Fig. 1.5 to the left). We have compared the acceleration 
predicted by this model with the acceleration observed in the experiment and 
find good agreement. 
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Fig. 1.5: Left: A toy model of the acceleration process. Electrons in the vicinity 
of the current channels are subject to an electromagnetic force, working 
to accelerate them along the ion flow. Crossing the centre of the channel 
the process reverses leading to an oscillating movement along the channel. 
Right: The normalised electron particle distribution function downstream 
of the shock. The dot-dashed line is a powerlaw fit to the non-thermal 
high energy tail. The inset shows a similar histogram for ion current 
density sampled in each grid point in the same slice as the electrons. 



To the right in Fig. 1.5 we have plotted the particle distribution func- 
tion for the electrons in a small slice in the box. We observe a powerlaw 
distribution. This should be understood as consequence of 1) the accelera- 
tion mechanism described above that directly relates the maximum kinetic 
energy of the electrons to the local ion current density and 2) the powerlaw 
distribution of the ion currents, as a consequence of the two-stream insta- 
bility, seen as an inset in the figure. The maximum acceleration observed is 
around f 7 ~ 80. Using the toy model and rescaling the ion to electron ratio 
of 16, used in the experiment, to the real value of 1836, we find the maximum 
energy gained by the electrons to be around 5GeV. 

The presented acceleration mechanism is essentially due to the electrons 
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oscillating in a potential, though as seen in Fig. 1.4 the true paths of the 
electrons are more complicated, and the radiative efficiency can be very high, 
because there are no free high energy electrons carrying away the kinetic 
energy such as in recursive acceleration scenarios. Moreover the properties 
of the process depend primarily on the local conditions of the plasma. 

In the chapter we estimate the thermalisation length for the ions, and 
find by extrapolating the fractional thermalisation observed at the boundary 
of the box, that the ions should thermalise in approximately 1500 ion skin 
depths. Using typical values for density in a gamma ray burst afterglow this is 
equivalent to 10^ m. We emphasise that the thermalisation length depends 
on the inflow velocity and mass ratio of ions to electrons among others, 
and a parameter study is necessary to uncover the true interdependence of 
parameters. 

Even though the two-streaming shock interface is estimated to be rela- 
tively thin, the high radiative efficiency implies that the non-thermal radia- 
tion observed in gamma ray burst afterglows and relativistic jets in general 
could be emitted from such a thin shell. 

1.5 The Global Structure of CoUisionless Shocks 

Collisions in "coUisionless shocks" are mediated by the collective electromag- 
netic field, and the scattering of the particles on the field slowly heats the 
particles. At some point the two-stream instability cannot be sustained, and 
the current channels become unfocused and decay, due to the thermal motion 
of the individual particles, which creates a warm turbulent medium with no 
significant large scale magnetic field. In Chapters 3 & 4 it is shown how 
magnetic field generation and particle acceleration are integral parts of rela- 
tivistic coUisionless shocks in the case of weak or absent large scale magnetic 
fields. 

To understand the impact on observations it is essential to investigate 
how far down stream of the initial shock that the two-stream unstable re- 
gion extends. With this in mind, in Chapter 5 I discuss the global struc- 
ture of coUisionless shocks. A range of experiments are presented, both 
three-dimensional models of pair plasmas and two-dimensional models of 
ion-electron plasmas. There is a fundamental difference between ion-electron 
shocks, where the mass difference leads to the ions dominating the dynamics 
and the electrons stabilising the ion channels, and a pair plasma, where the 
electrons and positrons form channels on the same timcscalc, and no shield- 
ing occurs. In the latter case the two-stream unstable region is significantly 
smaller than in the case of ion-electron shocks. 
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In the three-dimensional computer experiments we observe that the elec- 
trons and positrons thermalise fully, and the medium contains five different 
regions: The unperturbed upstream medium coming in from the left of the 
box; the first discontinuity in the velocity, with a two-stream unstable region; 
a warm thermalised region that is separated into a high and a low density 
state; another two-stream unstable discontinuity, where the warm shocked 
medium collides with the unperturbed downstream medium; and finally the 
unperturbed downstream medium. To verify that I have in fact resolved the 
full shock structure in a satisfactory manner, and the jump conditions have 
been established, I compare the experiment with a fluid simulation and find 
good agreement. From this experiment we can estimate that the two-stream 
unstable regions for electron-positron plasmas decay after 50-100 electron 
skin depths. 

In the second part of Chapter 5 I consider the global structure of ion- 
electron dominated collisionless shocks. With current computer capacities it 
is impossible to correctly model the global structure of an ion-electron shock 
in three dimensions. Two-dimensional collisionless shocks, being less costly 
computationally, remain a promising alternative, and I have investigated the 
applicability to understanding real three-dimensional models by performing 
large scale two-dimensional experiments (see Fig. 1.6), comparing them to 
the three-dimensional experiment discussed in Chapter 4. 

The particle distribution functions (PDFs) of the electrons for the two- 
dimensional and three-dimensional experiments are compared in Fig. 1.7. 
The slope indicated in Fig. 1.7 depends on the amount of heating in the 
upstream population, impacting the high energy part of the spectrum, and 
the down stream population, impacting the low energy part of the spec- 
trum. A warmer upstream population will be broader in phase space, and 
consequently the maximum is lower, giving rise to a steeper slope. The 
two-dimensional experiments have a slope index of 2.1, while the three- 
dimensional experiment has a slope index of 1.55. The difference in heating 
rates is understood in terms of the toy model, introduced above in section 
1.4 and discussed in Chapter 4, as a consequence of the different geometries. 

The physical significance of the the two-stream instability remains di- 
rectly related to the extent of the two-stream unstable region, and caution 
should be voiced about uncritically generalising results from two-dimensional 
experiments to three dimensions. My experiments seem to indicate that one 
will observe a faster thermalisation rate in two-dimensional experiment than 
what may be expected from three-dimensional experiments. 
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Fig. 1.6: The current density of the ions in a high resolution two-dimensional ex- 
periment. The dashed lines indicate the region used for constructing 
particle distribution functions. Length units are given in electron skin 
depths. 



1.6 A Next Generation PIC Code 

In Chapter 6 I present, together with C. Hededal, the first results from a new 
particle-in-cell code in development. The particle code that has been used 
to obtain the results described in Chapters 3-5 is limited to modelling the 
dynamics of charged particles under the influence of electromagnetic fields. 
In the new code, the concept of particles is generalised; most notably we 
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Fig. 1.7: Particle distribution function for the electrons in a slice indicated on 
Fig. 1.6. To the left is shown the PDF for the largest two-dimensional 
experiment, while to right the PDF for the three-dimensional experiment 
is shown. 



have introduced photons, and we consider microphysical interactions such as 
scatterings, decay, annihilation and pair production. 

Even though work still has to be done before we may start to investigate 
non trivial astrophysical scenarios, solid progress has already been made, and 
to test the infrastructure of the new code we have implemented Compton 
scattering as a simple scattering mechanism. The results are very promising; 
there is excellent agreement between theory and the numerical experiment. 

The new code will enable us to target problems that reside in the grey 
zone between the MHD and coUisionless plasma domains. This grey zone 
covers many astrophysical scenarios of great interest, among others internal 
shocks in gamma-ray bursts, solar flares and magnetic substorms, compact 
relativistic objects, and aspects of supernova remnants. 



2. GENERAL RELATIVISTIC 
MAGNETO-HYDRODYNAMICS 



Electromagnetic fields are ubiquitous ingredients in most astrophysical ob- 
jects. In the case of very compact objects or at cosmological scales, not only 
do electromagnetic fields interact with matter directly, but they also become 
a source of energy-momentum and impact on the metric curvature. Several 
general relativistic magnetohydrodynamics (GrMHD) computer codes have 
been developed and implemented recently for the study of compact relativis- 
tic objects and their surroundings [e.g. 3, 18, 19, 23, 40, 42], using both 
conserved and non-conserved formulations of the basic equations of motion. 
They are well-suited for their different purposes, but most of the implemen- 
tations above are designed for static space time backgrounds with diagonal 
spatial terms. 

In this chapter 1 present the analytic basis for and numerical implemen- 
tation of a code to solve the GrMHD equations. My approach is inspired by 
the pioneering work of Koide et al. [42] and related in spirit to the methods 
of Anton et al. [3] and Pons et al. [67]. Prom the beginning it has been de- 
signed to be general enough to solve the GrMHD matter evolution equations 
on any general time-dependent metric. This is an essential requirement if the 
code ultimately is to be coupled with numerical codes solving the Einstein 
equations, which evolve the metric. As far as the implementation is con- 
cerned I have currently implemented a fully parallelised 3D version of special 
relativistic MHD and a general relativistic extension of the hydrodynamics. 

In the following section 1 describe some of my motivations for developing 
the code. In section 2.2 I present the fundamental equations for GrMHD 
and adapt them to our specific approach. The equations are well known 
(e.g. [79]), but I make an effort to rewrite them in a form that is suited for my 
numerical purpose. Por clarity I first consider hydrodynamics and discuss the 
question of artificial viscosity and imperfect fiuids, to then extend the system 
to include electromagnetic fields. In section 2.3, I present the numerical 
algorithm that I have chosen to implement the equations with. Section 2.4 
contains a large test bed of demanding problems. Section 2.5 contains some 
astrophysics related tests of the code and finally, in section 2.6 I consider the 
crucial aspects of performance and scalability among others. 
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2.1 Motivation 

An important motivation for developing this kind of code is to make it pos- 
sible to study the evolution of cosmological magnetic fields in the primordial 
universe, taking into account the metric back reaction and coupling of grav- 
itational waves with the electro magnetic field. The WMAP satellite has 
already detected the first polarisation signal in the cosmic microwave back- 
ground radiation (CMBR) [63]. The Planck satellite and ground/balloon 
based experiments will improve the quality of the signal further in the com- 
ing years. Even though primordial magnetic fields make a very small con- 
tribution to the CMBR, in contrast to other imprints, they source vector 
perturbations and hence it may be possible to disentangle the weak signal 
from other sources through its unique character [29, 58, 66]. Turbulent pri- 
mordial magnetic fields can arise naturally during a phase transition, such as 
the transitions from an electroweak plasma and from the quark gluon phase 
to normal matter [80] . Alternatively, they may be produced during infiation 
[4] . If a signal from primordial magnetic fields is indeed detected, we would 
have yet another probe to understand early universe physics. Galaxies and 
clusters of galaxies at high redshift have been observed to contain magnetic 
fields comparable to present day galaxies. They have only rotated a few 
times during their short life, and this is difficult to explain without invok- 
ing primordial magnetic fields at some level. Dynamo theory alone does not 
seem to be enough [6, 29] . MHD simulations of turbulent helical fields have 
shown that an inverse cascade process operates which transfers small scale 
power to larger scales, changing the simple energy decay due to the expan- 
sion of the universe [15]. Until now, except from purely analytical analyses, 
the question of evolving magnetic fields in the early universe has primarily 
been tackled in two different ways. 1) Simple 3D turbulence experiments 
have been made, using existing non-relativistic MHD codes to address the 
possibility of inverse cascades which could alter significantly the longevity of 
large scale primordial fields; 2) Semi analytical arguments have been used 
to explore the couplings between primordial magnetic fields and the metric, 
neutrinos, effects from Silk-dampening, etc [46] If imprints in the cosmolog- 
ical microwave background from primordial magnetic fields are detected, it 
will be crucial to understand the evolution of the fields in a realistic manner, 
in order to constrain possible generation scenarios. 1 have verified the results 
by Christensson et al [15] using a purely special relativistic version of the 
code. With the code developed here, these questions may be addressed in 
a unified way, by performing large scale 3D experiments including general 
relativistic effects and couplings between the magnetic field and the metric 
perturbations. 
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Another strong motivation for developing a GrMHD code is the fact that 
it provides the perfect complement to the particle- and photon plasma codes, 
presented in the subsequent chapters, for the study of extreme astrophysics 
around compact objects and in jets. To understand the complex physics, 
we need to consider processes happening at many different time and length 
scales. A GrMHD code can be used to model the large scale dynamical 
flow and, as detailed in Chapter 6, provide realistic boundary conditions for 
microphysical studies of plasma instabilities and radiative processes. 

We note that the first results of couphng the full Einstein equations to 
the MHD equations has been published [20] only very recently, and that the 
field is still in its infancy. 



In numerical relativity it has proven very fruitful to exploit the so called 
3+1 split of the metric. Instead of working with a four dimensional man- 
ifold and the Einstein equations in the form of an elliptic nonlinear set of 
partial differential equations, an explicit split between temporal and spatial 
dimensions is imposed (though see [55] for an alternative four dimensional 
approach) . Assuming that we can construct a foliation of space time — usu- 
ally a very reasonable condition except maybe for near (naked) singularities 
— it is then possible to rewrite the Einstein equations as a hyperbolic set 
of evolution equations, some elliptic constraint equations and an associated 
Cauchy data set describing the initial conditions. This formulation lends 
itself easily to a numerical implementation and has been named the 3+1 
approach. 

The standard way of writing the metric in 3+1 form^ is: 



where a is called the lapse function, (3 is the shift or shear and 7 is the spatial 
3-metric. The contravariant version of the metric g'^'^ is written 



This form of the metric has the same number of degrees of freedom, namely 
ten, as in the obvious form g^i,. Here they are spread out as one for the lapse, 

^ Up to a plus or minus sign and a factor of for /3 
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2.2.1 3+1 Formulation of general relativity 



ds^ = -a^dt^ + {dx' + i3'dt) {dx^ + p^dt) , 



(2.1) 




(2.2) 
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three for the shear and finally six in the spatial curvature. Therefore, any 
metric which is not null may be written in this form. 

In this thesis I only consider the evolution of matter and fields in a back- 
ground space time although through the Einstein equation they are sources 
for the metric fields. Thus, it is important to leave a, j3 and 7^, unspecified, 
making the design ready for integration with evolving metric fields. 



2.2.2 Different coordinate systems 

The global coordinate system Eq. (2.1) is often called the star fixed coordinate 
system (SFCS), because in most applications it is asymptotically fiat and, 
therefore, connected to inertial observers at infinity. If we consider instead 
local observers who do not observe any shear and measure time in terms of 
local clocks, their line element must be given as 

ds'^ = -dp + -fijdx'dx^. (2.3) 

This coordinate system is denoted the local laboratory frame (LOLA frame), 
and I write any quantity in this coordinate system with a hat. In the LOLA 
frame in many interesting cases 7^^ is almost diagonal and one could then 
easily rescale the problem as done by Koide et al. [42] to evolve matter and 
fields as seen by local observers, or FIDOs^ instead. 

I have done so but to keep my approach general, I have exploited the 
idea to always rescale the diagonal in the metric, even though it may well 
be non diagonal. Because the off diagonal terms in the spatial part of the 
metric often are comparable in size to the diagonal ones, I have effectively 
normalised the metric. Since the metric is almost a FIDO metric I have 
named it the pseudo FIDO frame (PFIDO) frame. In this frame the metric 
tensor is given as 

ds'^ = -dp + %jdx'dff , (2.4) 

and there are only three non-trivial terms in the PFIDO metric, because all 
but the non diagonal terms in Eq. (2.4) have been normalised. 

The central idea of my numerical scheme is to use the PFIDO frame to 
measure all physical quantities. The PFIDO frame is only defined locally and 
we still need to use the global coordinates connected to the SFCS to measure 
distances. The general way to construct an equation is first to derive it in the 

^ FIDOs are fiducial observers whose metric is defined as that seen by observers in local 
inertial frames. 
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SFCS and then to transform the tensors and vectors from the SFCS to the 
PFIDO frame, while keeping the derivatives and differentials with respect 
to the SFCS. It is central that the transformation from the SFCS to the 
PFIDO frame is completely linear and simple, even for generally evolving 
coordinates. Had we, instead, chosen to go all the way to a FIDO frame in 
the general case, we would have had to invert a matrix, the metric, at every 
point for every time step. The PFIDO frame is a healthy in-between, which 
gives us almost all of the advantages of the FIDO frame but at a much lower 
cost. 

Intuitively it is clear that when going to a local frame of reference the 
curvature of space only manifests itself as extra external Coriolis-like forces, 
giving some extra terms in the evolution equations below. From a numerical 
view point there is an added benefit when we consider space times with a 
strong shear or frame dragging, ie. points where /3' is large. The standard 
example of this is the Kerr metric in Boyer-Lindquist coordinates. Inside the 
ergosphere, from the point of view of the SFCS, everything is rotating around 
the black hole in the same direction as the spin of the hole. The closer we 
are to the event horizon, the faster the rotation induced by the shear. From 
a local observers point of view in the PFIDO frame though, there is no shear 
and the locally defined velocity is much smaller. The locally defined velocity 
is the truly interesting velocity, since it arises due to physical processes, while 
the apparent high velocity seen by an observer attached to the SFCS is partly 
due to the geometrical structure of the background space time and partly due 
to physical processes, thus, a result of the chosen reference frame. Near the 
horizon the shear-induced frame dragging velocity can be much greater than 
the local velocity, and we can run into problems with numerical cancellations 
smearing out variations in the local velocity. Yet this is avoided by choosing 
to work in the PFIDO frame. 

From the hue elements Eq. (2.1) and Eq. (2.4) we may derive the trans- 
formation laws. In particular we have 

adt^di, (2.6) 
^/^i{dx' + (3'dt) = dx' , (2.7) 

and coordinate differentials are contravariant vectors. Then, any contravari- 
ant vector Lf^ transforms like 

= aU^ , U' = V7i7 {U' + . (2.8) 

It is a matter of linear algebra to show that covariant vectors transform like 

Ut = -{Ut-P''Ui) , U, = -^U,. (2.9) 
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Tensors transform as the product of vectors by their very definition. We refer 
the reader to App. B for a complete hst of transformation relations that have 
proven useful when deriving the equations in this chapter. 

2.2.3 Basic equations 

The basic fluid equations follow from conservation laws. The conservation of 
the baryon current gives 

V^(pC/'^) = 0, (2.10) 

where is the covariant derivative, p is the rest mass density and If^ is 
the four velocity in the SFCS coordinate system. The conservation of the 
energy-momentum tensor leads to a similar expression 

V^TIf = 0. (2.11) 

The version that we have chosen to use of the energy-momentum tensor for 
a fluid is given as 

Tl^HD)^ = PhU'U, + - , (2.12) 

where P is the pressure, h—l + eint + P/p is the relativistic enthalpy, eint is 
the internal energy, and rjaf^ is the shear viscosity. It has the deflnition 

^ 1 (h^^S/^U'' + W^Vo^Ui") , (2.13) 

where h^'^ projects into the fluid rest frame 

hi"^ ^ Ui'U'' ^ gi'r (2.14) 

We consider the energy-momentum tensor in mixed form as the basic hy- 
drodynamical object to evolve, because even for general metrics the pressure 
term disappears in Eq. (2.12) for off-diagonal components [27]. This is not 
the case for the purely co- or contravariant versions. 

The energy momentum tensor of the electromagnetic fleld is 

^(W). = - \KF-''F,a , (2.15) 

where F'^^ is the electromagnetic fleld strength tensor. 



2.2 The GrMHD equations 



19 



We can simplify the covariant derivatives significantly by using the fol- 
lowing identities 

VJU^ = -^^=^9, (yqp/C/M) , (2.16a) 
V^TT = (v^HM^.^) - iT^'^d.g^ , (2.16b) 

= ^7^^^^ (a/HM^/) , (2.16c) 

where / is a scalar function, a vector, T^^ is any symmetric tensor, F^^'^ 
any antisymmetric tensor and is the determinant of the metric. 

2.2.4 Selecting evolution variables 

We have chosen our field variables with respect to the PFIDO frame and the 
basic evolution variables take the form 

D ^-fpU^ ^-fpW , (2.17) 
E = -7ff - D = 7 {phW^ -P- pW) , (2.18) 
Vi = ^/lUlfl = ^/lUlphWUi , (2.19) 

where W = is the Lorentz factor of the fluid with respect to the PFIDO 
frame and 7 = \/\\'^\ \ is the square root of the determinant of the spatial 
metric. Looking at Eq. (2.16) and Eq. (2.8) we see that the reason for 
choosing the factor 7 in front of the variables in Eqs. (2.17)-(2.19) is to 
cancel out a/— | If?! | in Eq. (2.16). The subtraction of the relativistic mass 
density in the definition of the total fiuid energy density is done in order 
to cancel the rest mass energy density, which could otherwise jeopardise a 
numerical implementation when the flow is non-relativistic. 



2.2.5 Hydrodynamic equations 

In order to highlight the physical content I flrst write down the equations of 
motion in the case where there are no electromagnetic fields: Ti^" — T^^j^y 
To find the equations of motion, we use Eqs. (2.8)-(2.9) and their extension 
to mixed typed two-tensors (see App. B) together with the rules for covariant 
derivatives Eq. (2.16) and the fundamental equations of motion in the SFCS 
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Eq. (2.10) and Eq. (2.12) 

dtD = -djDv^ , 
dt [S + S*] 



-a 



1 

+ - 

a 



{S + 7P) + Sw, 

- [DhW {dt + v^dj) + 7P {dt - P^dj) 

+Eldt + ^ldj In a 

- d, (7P/3O + {P'M. - Mt) , 



dt [Pi + E*] = -dj 



- di [a'-fP] + aMi , 



where the normal three-velocity has the usual definition 



1p 



the transport velocity is the three velocity seen from the SFCS 



a 



P JJt ' 



the geometrical terms are 
and the viscosity terms are 



1 





-1^1 , 
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n = 
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HVt = 
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a 



(2.20) 



(2.21) 
(2.22) 



;2.23) 



(2.24) 



(2.25) 

(2.26) 
(2.27) 

(2.28) 
(2.29) 



Even though the evolution equation for the energy has become a bit more 
comphcated than in the special relativistic case (a = 7 = ^/7i7 = 1, = 0), it 
represents a substantial simplification in that relations between the different 
variables reduce almost to the special relativistic form. Hence for example 
the Lorentz factor W may be computed as W = [1 + ^ijWU^Y^'^ bearing in 
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mind that the diagonal is already normalised. Let us consider a space time 
without any off-diagonal spatial components but with an arbitrary shear. For 
example Boyer Lindquist coordinates in an extreme astrophysics context or 
the uniform curvature gauge in a cosmological context. In these examples, the 
shear viscosity is identical to the special relativistic form. This is because the 
PFIDO frame reduces to a FIDO frame of reference. To handle coordinate 
systems that penetrate the event horizon of a black hole, for example the 
Kerr-Schild coordinates, we need at least one off-diagonal spatial component 
[17]. In this case extra terms in the shear tensor arise, but changes are 
minimal. 



2.2.6 Artificial viscosity 

It was argued by Anninos & Fragile [2] that in order to make a consistent 
relativistic finite difference code with artificial viscosity (AV) it is crucial to 
use a viscosity that has been defined in a physically sensible manner, oth- 
erwise it will break down for flows with high Lorentz factors. An efficient 
AV should be covariant in its deflnition, such that, the code can easily be 
adapted to general relativity, be physically meaningful, respect energy con- 
servation, and reduce to some normal Newtonian AV formulation in the non 
relativistic limit. We know of no implementation so far that has respected all 
of the above points. Indeed it seems that the prevalent thing is constructing 
a mock-up "viscous pressure" using the prescription P — > P -|- Qvisc and then 
include a directional dependence such that the effective energy-momentum 
tensor takes the form 

Tj^;;^^ = {ph + Q^iscWU^ + g^^P + Q^"" . (2.30) 

Such a viscosity may be able to deal with mildly relativistic shocks but it 
does not even reduce properly in the non relativistic limit. 

A general imperfect fluid energy-momentum tensor may be written 

Tl^HD) = phUW + g'^^'P + Q"" , (2.31) 
Q^"" = -2ria^"' - iOh^" , (2.32) 

where r] and ^ is the shear and bulk viscosity coefficients, 9 = V^lf^ is the 
expansion of fiuid world lines, and cr^'^ is the spatial shear tensor (see Eq. 
(2.13)). In the non relativistic limit we find that 

T''-D^ ^ p^,^^ ^ ^2.33) 
T'' pv' , (2.34) 
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which shows that any consistent shear viscosity should reduce as 

{Q'\Q'') ^ {v'nj,nj) (2.35) 
Tij = lyij {div^ + djv') (2.36) 

in the non relativistic hmit. Here Uij is some viscous constant, which could 
depend on the numerical grid spacing dx\ the local sound speed and other 
factors. Neither the viscous pressure formulation (Eq. (2.30)) nor the bulk 
viscosity ^Oh^'' reduce properly in the limit. Only the shear viscosity -qa^" 
does so. The shear viscosity is included directly in the energy-momentum 
tensor and it is by construction covariant and preserves energy and momen- 
tum. 



2.2.7 Electromagnetic fields 

The 3-1-1 formulation of Maxwell's equations was originally calculated by 
Thorne & MacDonald [79] and may be written (see also Baumgarte & Shapiro 
[8]) 

dnE' = 47r'ype, (2.37) 
daB' = , (2.38) 
daE' = e'^^dj{a^Bk) - Aira^f + dj [jS^^E' - (3'-fE^] , (2.39) 
daB' = -e'^'djia'^Ek) + dj [(3^'yB' - P'-fB^] , (2.40) 

where E^, B^, pe and J* are the electric field, magnetic field, charge density 
and current density as seen by observers in the SFCS frame. With the goal 
of simplifying the equations, we absorb the determinant of the 3-metric in 
the definition of the different fields. Furthermore we use the fields as seen by 
observers in the PFIDO frame. The Maxwell equations then become 



9.r = 47rp,, (2.41) 

diB' = , (2.42) 

dtS' = e'^^dj{aBk) - 47ra7 + dj [(3^£' - (5'£^] , (2.43) 

= -e'^^dj{a£k) + dj [f3^B' - (3'B^] , (2.44) 

where B' = -^B' = jB\ Bi = J^aBi = jBi, £' = -^E' = jE\ 

^ / Yl 1 V ^ / It 1 



£i = \/li'ilEi = '~fEi, pg = and J = -^=J^. Except for the shift terms 
and some lapse factors, this equation set is identical to the special relativistic 
Maxwell equations. 
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The energy and momentum equations are modified in the presence of 
electromagnetic fields, refiecting the transfer between fields and fiuids. 



where is the four current vector. After some algebra we find 
dtE^ 



(2.45) 



+ 7 WF^,J^ - Ft^J^] = . . . + ^rSi 



a- 



dtV^ 



7 



. . . + 



a 



7 L 



(■ijkJ'' B^ 



+ 



a 



7 



JxB + p^-£ 



(2.46) 



(2.47) 



It is worth noticing that the result practically reduces to special relativity 
except for the pref actor a^~^. 

2.2.8 Ohm's Law 

If we consider relativistic MHD, we have to supply an Ohm's law to link the 
electric and magnetic fields with the current density. A relativistic version 
of the standard non-relativistic Ohm's law may be written [8, 47, 54] 



where rjc is the resistivity. Using Eq. (2.8) it reduces to 



(2.48) 



W ' 
w ' 



hi 
1 



hi 



(2.49) 



Except for the Lorentz factor W and the single geometric factor, this is 
identical to the standard non relativistic result. 

In this thesis the ideal MHD condition will not be used directly, since 
resistivity is applied in the code. However, taking rj^ = and assuming the 
ideal MHD condition Faradays law Eq. (2.44) in the SFCS may be reduced 
to [8] 

dt-iB' = dj {{UyUhB^ - {UyW-fB') , (2.50) 
which in our notation is 



dtB' = dj {jfB^ - v^B') . 



(2.51) 
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2.3 The Numerical Algorithm 

I have used the equations of motion Eqs. (2.20), (2.46), (2.47), (2.44) together 
with Eq. (2.43) for the current density and an Ohms law Eq. (2.49) as a 
basis for the general relativistic code, but even though many mathematically 
equivalent forms of the equations of motion exist, they may lead to numerical 
implementations with radically different success rates. In this section, I detail 
some of the concepts I have used to deal with the problems that inevitably 
arise when solving a set of equations numerically. 

The most important choice is to determine if we want to exploit the 
characteristic structure of the equations or just directly use finite differencing 
to solve the equations. In keeping with the tradition in Copenhagen I have 
chosen the latter. This has helped to develop the code in a relatively short 
time span and I am indebted in my reuse of techniques and tricks from the 
non relativistic codes developed in Copenhagen. 

The next fundamental choice is the form of the equations. Either we can 
use a flux conservative or a non conservative formulation. There are benefits 
to both: In the fiux conservative formulation, the Rankine-Hugoniot jump 
conditions are automatically satisfied across shock fronts even if the model 
does not resolve the shocks entirely. This is not the case for a non conser- 
vative formulation. On the other hand: In a fiux conservative formulation, 
one of the conserved variables is the total energy. It contains contributions 
both from the fluid and from the electromagnetic flelds. If the plasma is 
strongly dominated by the electromagnetic fields, the internal energy, the 
difference between the total and electromagnetic energies, can be swamped 
by numerical noise and round off. Another problem — albeit technical — 
is that the conservative variables in the MHD case are related algebraically 
to the so called primitive variables through a sixth order polynomial. There 
is no analytical solution to the problem, and an expensive numerical root 
finder method has to be used. 

I have chosen a cross breed solution: 1 use conservative variables for the 
hydrodynamics, while in the case of MHD, I do not include the magnetic 
energy and momentum in the total energy £ and covariant momentum Vi- 
The basic reason for not using conservative variables is due to the problems 
with magnetically dominated plasmas. As an added benefit, 1 circumvent 
the problems of finding primitive variables through non analytical methods. 
Nonetheless, still at every time step it is necessary to find the four velocity 
tj^ and enthalpy h from the total hydrodynamic energy S and covariant 
momentum Vi- 
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2.3.1 Primitive variables 

Given the dynamical variables D, £ and Vi in Eqs. (2.17)-(2.19) together 
with the equation of state for an ideal gas 

P = (r - l)pe,„i = ^^p{h - 1) , (2.52) 

where F is the adiabatic index, I define two derived quantities 

X^^^{h-l)W + W-l-^^^, (2.53) 

Y='^^h\W^-l). (2.54) 

Using Eq. (2.54) to solve for W and inserting the solution into Eq. (2.53) 
a fourth order polynomial va. hm = h — 1 may be constructed, which only 
contains X, Y and F in the coefficients, viz. 

hi + 2[F + 1] Z,^ + [1 + r(4 - FX(2 + X) + 2Y)] hl+ 

[1 - TX{X + 2) + (1 + T)Y] h„,+ 

r^{l + Y)(Y -X"^ -2X) = 0. (2.55) 

When the desired root has been found, it is trivial from Eq. (2.54) to obtain 
WUi = — 1 and then any other desired quantity. Fourth order polynomials 
may be solved iteratively using a range of different root finder methods, such 
as the Newton-Raphson method. I tried this, and even though it most often 
worked fiawlessly and was fast, for certain corner cases, it is both unstable 
and slow. Slowness in a few cases may be acceptable, but if the method 
crashes, the simulation crashes. Stability is the key. An alternative is to use 
an analytic formula for the roots, but great care has to be taken. In any naive 
implementation, for example taking directly the output from Mathematica, 
the coefficients will cancel numerically at the slightest difference in scale of 
the four velocity and the Lorcntz boost and the result will be imprecise. In 
the end I settled on a method detailed in [1] to reformulate the problem in 
terms of roots in one third order and four second order polynomials. I find the 
roots using stable formulae, which guard for cancellations, from [69]. With 
this approach the code runs most tests using single precision variables and 
only for the most extreme cases (high Lorcntz boost and very low pressure), 
we have to fall back to double precision. The solver is not only rock solid 
but also very fast. Properly implemented with no if-lines and all calculations 
vectorised, it takes approximately 20% of a time step, and therefore does not, 
in any way, dominate the problem. Note that a related approach has been 
reported in [19]. 
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2.3.2 Artificial viscosity 

I do not try to solve, neither exactly nor approximately, the Riemann problem 
at cell boundaries. Instead, I use finite difference derivatives. To stabilise the 
algorithm it is critical to add AV. During the development of the code I have 
tried many different formulations both inspired by the non relativistic codes 
developed in Copenhagen, classical formulations of AV and the self consistent 
AV detailed in [2] . In the end I settled for an AV based on a physical model 
of shear viscosity derived from the energy momentum tensor of an imperfect 
ffuid (see section 2.2.6). To determine the viscosity coefficient r] in front of the 
shear viscosity in Eq. (2.12) I have extended the prescription already used in 
the non relativistic codes in Copenhagen [61], and use a Richtmeyer-Morton 
type hyper viscosity that depends on the local conditions in the fluid: 

rjij = Axij viCs + i'2,\v\ + P2Al\dJj^\<o , (2.56) 
I^Xij = ]^Dh [Ax^ + Ax^] , (2.57) 

where Cs is the relativistic sound speed. A/ = max(Ax*) and | ■ |<o means that 
the strong shock viscosity only is operative where there is a compression of 
the fluid. Except for the sound speed, the only other changes in the coefficient 
compared to [61] are the use of Dh, as seen by an observer in the local 
PFIDO, frame instead of the mass density p, and the use of the divergence of 
the four velocity in the relativistic case compared to the normal divergence 
of the spatial three velocity in the non relativistic case. It is non trivial to 
find the time derivative of the Lorentz boost W . We found by experimenting 
with different, mathematically equivalent prescriptions, that by far the most 
stable formulation is 

dtW = ^dtU% . (2.58) 

The shear viscosity, given in Eq. (2.13), contains time derivatives of the four 
velocity too. In the code I use a third order Runge-Kutta integrator for 
the normal dynamical variables. I evaluate the four velocity derivatives by 
explicit derivatives, storing old velocities three sub time steps back in time. 
This way I get third order correct time derivatives. Unfortunately they are 
not correctly time centred and I speculate that some of the problems I see 
in the test problems below for high Lorentz boosts may be due to the time 
derivatives lagging approximately half a full time step compared to the rest 
of the terms. In the energy and the momentum equations (2.21) and (2.22) 
AV terms arise both on the right hand side and in the time derivative. I have 
currently not included the time derivative of the shear viscosity in the code. 
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2.3.3 The magnetic field 



The equations are evolved on a staggered mesh (see below) and the divergence 
free condition Eq. (2.38) of the magnetic field is naturally conserved. To raise 
the entropy in magnetically driven shocks I use the exact same formulation 
as in [61] for the resistivity rjc since the Maxwell equations by construction 
comply with special relativity, and the only change has been to substitute a 
relativistic correct expression for the fast mode speed. 

Ohms law Eq. (2.49) and Amperes law Eq. (2.43) are used to derive the 
electric field and the current density respectively. We use an explicit time 
derivative to evaluate the displacement current. Even though it is lagging 
behind with half a time step, like the time derivatives of the four velocity, 
it has proven very effective in limiting the magnetically driven wave speeds 
except when the Alfven velocity becomes close to the speed of light. The 
magnetic part of the code is calculated following the scheme 

• Calculate the resistivity rjc- It is proportional to ubJ^s- 

• Estimate the electric field: £* = -^=v^ x 



• Calculate S*^ and find the displacement current using an explicit time 
derivative. 

• Calculate an estimate for the current aJ*^ — e^^^dj{aBk)-\-dj \J3^S*^ — P^S*^] . 

• Lower the current and find the final electric field £i = + £*. 



• Proceed calculating Faradays law Eq. (2.44) and the energy and mo- 
mentum sources Eqs. (2.46) and (2.47). 

I have tested different variations of the scheme above using the full version of 
the current density j\ including the displacement current, to find the final 
electric field. Even though formally better, it turned out to be less stable, 
giving short wave oscillations and essentially the same results. 



I have implemented an MHD version of the above equations, currently re- 
stricted to special relativity. A pure HD version has been made for general 
relativity with diagonal metrics. To test the code, I have applied a battery of 
tests that are presented below. In all tests I have used a 3 dimensional version 



Use the displacement current to update the current J — J 
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of the code. The boundary conditions are implemented in the y direction by 
design, and therefore our ID box has the size (1, Ny, 1). If not stated other- 
wise, in all runs, the weak shock viscosity coefficients are i/i — i/^ — 0.029, 
the strong shock viscosity coefficient is 1/2 — 0.55, the magnetic resistivity 
coefficient (see [61]) is z/^ = 1 and the Courant limit is Cat = 0.3. The code 
can handle more extreme problems by tuning the different numbers, but I feel 
that it is important that the code "just works" ; in real physical applications 
the results should not rely too much on the tuning of these technical param- 
eters, since that would question the validity of the results. As an example, 
by just decreasing the Courant limit and the weak viscosity z/i I am able to 
run the wall shock test with a Wmfiow = 5 and obtain satisfactory results. 
Only in two of the magnetic tests, I have tuned the coefficients to facilitate 
the comparison with other codes. 

2.4.1 Hydrodynamical tests 

The code has been developed without extending any preexisting relativistic 
fluid dynamics code, and it is important to demonstrate that it can solve 

correctly a variety of purely hydrodynamical problems. Fortunately, the 
analytic solution to hydrodynamic shock tubes is known [49, 68, 78]. I have 
used the RIEMANN program published by Marti and Miiller [50] to generate 
the analytic solutions. 

Blast waves 

The blast wave is a problem with two domains initially at rest with a discon- 
tinuous jump in the density and pressure. A blast wave is launched at the 
interface with a very thin shell of high density. The fluid separates in five 
different states. Two initial states at the left and right boundary, a rarefac- 
tion wave, the contact discontinuity and a shock wave. This setup is ideal 
for testing how diffusive the scheme is, since the shock wave, for suitable pa- 
rameters, is very thin. The initial states for the three problems we consider 
are shown in Table I. 



Table I 


Problem I 


Problem II 




Problem III 


Blast waves 


Left Right 


Left Right 


Left 


Center Right 


Pressure 


13.33 0.001 


100 0.1 


100 


0.01 100 


Density 


10 1 


1 1 


1 


1 1 


Gas Gamma 


5/3 


1.4 




5/3 



Problem I, shown in Fig. 2.1, is a classic shock tube, that most relativistic 
codes have been tested against (see [50] for a compilation). Ideally the right 
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Fig. 2.1: Problem I: A mildly relativistic blast wave problem. Notice the slight 
oscillation at the edge of the shock front. This is due to the large jump 
in pressure at that point. 



state should have zero pressure but due to numerical reasons we have set it 
to 0.001. A small weakness of the code is already visible in this test. When a 
high mass density region separates from a low density region, such as at the 
contact discontinuity in Fig. 2.1, there is a certain amount of stickiness. It is 
in fact a feature to avoid low density regions to develop into true vacuums, 
making the code crash, but it also makes advecting high density blobs become 
more diffusive at the trailing edge. The shock velocity is maintained to a 
very high precision and the rarefaction wave is near perfect too, even at low 
resolutions. 

Problem II is a more relativistic variation of problem I. The shock wave 
is propagating with 0.92c. At t = 0.4 the shell has a thickness of Ay = 0.023 
or 11 grid zones at a resolution of 500 points. The AV spreads out the 
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Fig. 2.2: Problem II: A relativistic blast wave problem. Our code has some prob- 
lems with maintaining sharp contact discontinuities at points, where a 
high density blob is moving away from a low density area, such as just 
behind the high density shell. 



discontinuity over 6 points and this explains why the shock wave is under 
resolved at this resolution. 2000 points are needed to get a reasonable solution 
at t = 0.4. Notice also, that the diffusion in the density impacts on the flat 
profiles of pressure and velocity. 

Problem III is the most extreme shock tube. To make a different setup I 
have removed the rigid boundaries and instead imposed periodic boundaries 
(see Fig. 2.4). A similar problem was considered by Marti and Miiller [49]. 
Compared to problem II, the pressure in the right zone is also lowered, and 
the equation of state is more sensitive to the pressure. 

When the two shock waves collide ait = 0.26 a very dense shell is created. 
To track the evolution in an easy way, I have plotted the maximum density 
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Fig. 2.3: Problem III: Colliding blast waves. The evolution of the maximum in 
density as a function of time is shown. A resolution of 8000 points is 
needed to resolve the very thin shell of high density that is created when 
the two blast waves collide, and to accurately calculate the post shock 
profile, while with 2000 points we marginally resolve the preshock solution 
at t = 0.26. 




Fig. 2.4: Problem III: The solution at t = 0.2, just before the two shock waves 
collide. 
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Fig. 2.5: Problem III: The system, at the colhsion at t = 0.265. Notice we have 
changed the scale of both the x- and y-axis to reflect the the large change 
in density, and visualise the thin structures. See Fig. 2.4 for legend. 




Fig. 2.6: Problem III: The system, after the collision at t = 0.3. See Fig. 2.4 for 
legend. 
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Fig. 2. 7: Problem IV: The wall shock problem. The solution is shown at t = 2 and 
the resolution is 200 points. 
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Fig. 2.8: Problem IV: The same as in Fig. 2.7, but the resolution is 400 points. 

Notice that the number of points in the shock interface stays the same 
for different resolutions; about 3, 2 and 11/2 points for the different 
velocities. 



as a function of time in Fig. 2.3 for different resolutions. To resolve the 
preshock state reasonably well, at least 2000 points are needed, while 8000 
points are necessary to resolve the high density region and the post shocks. 
In [49] 4000 points were needed using a shock capturing PPM method to 
accurately model their problem. 



The wall shock 

The last hydrodynamical problem I have tested against is the wall shock. A 
cold fluid comes in from the right and hits a wall at the left edge where it 
is reflected. The inflow density is p = 1 and the adiabatic index is F = 4/3. 
When reflected a warm dense medium builds up. Figs. 2.7 and 2.8 shows 
the solution at different resolutions and time t = 2 for mildly relativistic ve- 
locities of Vs = 0.9 and downwards. The analytic solution to the wall shock 
problem may be found in [2] and [50]. 



It is clear from the above tests that the code is working very well up to 
a Lorentz factor of about W = 2.5. For higher Lorentz factors the current 
artificial viscosity implementation becomes problematic. I believe there are 
two problems with the current implementation: We use explicit time deriva- 
tives for the four velocities, but exactly because they are explicit, for a given 
time step t they are found at t — ^dt, and if the fluid is highly relativistic this 
will make a difference. In the wall shock, I observe that only decreasing the 
Courant limiter from the stock 0.3 to 0.01, I can reach an inflow velocity with 
a Lorentz boost of 3.5. Anninos & Fragile [2] have developed, to our best 
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knowledge, the only explicit AV based code that can handle high Lorentz 
factors. This is possible, because they include the time derivatives of the 
viscosity. 

2.4.2 Magnetohydrodynamical tests 

To validate the magnetic aspects of the code, I have performed a range of 
tests. Unfortunately, in relativistic MHD, no analytic solution is known to 
the Riemann problem, and I have to rely on comparison with tests considered 
by other groups using different codes and methods. Komissarov published in 
1999 a testbed [43] (hereafter K99) with different shock tubes. Unfortunately 
there were some errors in the tables, which are corrected in [44]. Some of 
the tests were used by De Villiers and Hawley [18] and Gammie et al [27] 
to validate their respective GrMHD codes. I have continued this trend by 
performing the same tests as in [18]. They augmented the testbed of K99 
with an Alfven pulse test that tests for correct wave speed of Alfven waves 
at different background fluid speeds and degrees of magnetisation, and a 
more complete set of magnetosonic shocks. Presented below are tests of 
magnetosonic shocks, magnetised shock tubes and similar Alfven pulses. 

Magnetosonic slioclis 

In Fig. 2.9 I present a collection of four different standing magnetosonic shock 
waves. The parameters of the different waves may be found in table II and 
have been taken from [18]. In the most extreme shock, the Fast Shock III, 
we had to decrease the Courant limit to Cat — 0.1 and the shock viscosities 
to {ui,U2) = (0.001,0.03). 

For all cases the solution is in excellent agreement with the analytical 
solution. Only in the case of the slow shock a slight over density has built 
up and is propagating away from the shock wave. This might be due to a 
relaxation of slightly imperfect initial conditions, and the solution has instead 
settled to a new static solution with a small difference in the parameters. In 
the cases of the fast shocks initially there is a perturbation too, but only as a 
small temporary ripple. In the cases of the Fast Shock II and III (the lower 
plots in Fig. 2.9) the ripple has already been advected out of the box, while 
in the case of the Fast Shock I it can still be seen at the right edge of the 
figure. 

Magnetised shock tubes 

I have performed two magnetised shock tube tests and the parameters may 
be found in table II. The first is a relativistic version of the classic shock 
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Fig. 2.9: Problem V: Magnetosonic shocks. The slow shock is top left, fast shock 
I is top right, fast shock II bottom left and fast shock III is bottom right. 
The buildup in the right side of the slow shock is due to interaction with 
the boundary. In the other shocks, the solution is close to perfect and 
buildup does not occur. 
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Fig. 2.10: Problem VI: Magnetised shock tubes. To the left is the relativistic 
version of the Brio & Wu shock tube, to the right the K99 shock tube 2. 
Compared to Figs. 6 and 7 in [18] and Fig. 6 in [43] it is clear that most 
of the different waves have the correct amplitude, but there are problems 
with too high wave speed and therefore errors in the rarefaction wave. 
This is most pronounced for the K99 shock tube to the right. 
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Fig. 2.11: Problem VII: Alfven pulse test. We start two Alfven pulses at y = 1.5. 

The wave speeds depend on the background fluid velocity and the degree 
of magnetisation. We begin to get significant errors when va ^ 0.7c. In 
all figures, the time is selected to have the two waves line up. This is 
not the case for ALFl and ALF3. 
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Table II 


Slow Shock {Vs = 0.5) 


Fast Shock I {Vs = 0) 




Left 


Right 


Left Right 


Pressure 


10 


55.33 


2.015 5.135 


Density 


1 


3.322 


1.406 2.714 


Four Vel 


(0,1.53,0) 


(0,0.957,-0.682) (0,1.78,0.114) (0,0.922,0.403) 


Mag Field 


(0,10,18.28) 


(0,10,14.49) 


(0,3.33,2.5) (0,3.33,4.52) 


Gamma 






4/3 


^final 




2.0 


2.5 


Grid size 




512 


1024 




Fast Shock II {Vs = 0.2) 


Fast Shock III {Vs = 0.2) 




Left 


Right 


Left Right 


Pressure 


2.015 


2.655 


2.015 34.99 


Density 


1.406 


1.725 


1.406 8.742 


Four Vel 


(0,1.78,0.114) 


(0,1.479,0.28) 


(0,3.649,0.114) (0,0.715,0.231) 


Mag Field 


(0,3.33,2.5) 


(0,3.33,3.25) 


(0,3.33,2.5) (0,3.33,6.52) 


Gamma 






4/3 


^final 




2.5 


2.5 


Grid tiizc 




1021 


512 




Relativist 


;ic Brio & Wu 


Shock tube 2 from K99 




Left 


Right 


Left Right 


Pressure 


1.0 


0.1 


30 1.0 


Density 


1.0 


0.13 


1.0 0.1 


Mag Field 


(0,0.75,1.0) 


(0,0.75,-1.0) 


(0,0,20) (0,0,0) 


Gamma 






4/3 


^final 






1.0 


Grid size 




2048 


512 









Tabic III: Alfven pulse tests 








Test 


(3 




v+ 


V- 10"^ X A+ 


10^ X A- 


time 




ALFl 


0.05 


0.0 


1.04(0.85) 


-1.04(-0.85) 5.0(5.0) 


5.0(5.0) 


1.17 


7.8 


ALF2 


0.315 


0.249 


0.48(0.47) 


0.00(0.00) 4.8(4.7) 


5.2(5.3) 


2.13 


5.4 


ALF3 


0.1 


0.8 


1.09(0.95) 


0.33(0.34) 3.8(2.5) 


6.2(7.5) 


1.46 


10.7 


ALF4 


0.315 


0.088 


0.33(0.33) 


-0.17(0.17) 0.49(0.49) 


0.51(0.51) 


4.04 


5.0 
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tube of Brio and Wu [12]. The shock tube is not very extreme and with a 
resolution of 2048 points, just like in [18] we clearly resolve all shock waves. 
The solution is shown att — 1. Comparing with Fig. 7 in [18] we see that the 
wave speeds are wrong. The right rarefaction wave has reached y — 1.1 and 
is superluminal while it should have reached y = 0.9. The left rarefaction 
wave is in good agreement with Fig. 7 in [18] propagating with v ^ 0.68. 

1 was only able to obtain a stable solution of shock tube 2 of K99 by 
lowering the viscosity to i/i — 0.001 and enhancing the magnetic resistivity 
to ub — 3.0. The shock tube is on the limit of the codes capability; small 
oscillations in the density p just behind the shock wave in Fig. 2.10 are 
evident and there are large errors in the rarefaction wave, which propagates 
superluminally at f = 1.2c. The forward going shock wave is only slightly 
wrong propagating with v ^ 1, where it should be going with v — 0.95. 



Alfven Pulse test 

The test is conceptually very simple. In a background with constant magnetic 
field and velocity in the y direction we set up a small square pulse in the 
perpendicular velocity component v^. It splits into two waves that travel 
with the Alfven velocity. The test is presented in [18] and although simple 
in concept it will easily reveal any errors in the wave speed. Since we use a 
direct finite difference technique to solve for the magnetic field it is critical 
that the displacement current is calculated correctly when the Alfven speed 
approaches the speed of light. It is already evident from the shock tube 
test above that this is not always the case, and this test has been invaluable 
during the implementation, for assessing different schemes to calculate the 
displacement current. 

Initially there is a constant background magnetic field and a constant 
background fluid velocity v^. On top of that a small square pulse with 
transverse velocity is superimposed. The pulse will split in two waves 
travelling with the Alfven velocity, given by [18] 



(,59) 



where = |6|^/(p/iW^^) and is the magnetic fleld measured in the fluid 
rest frame. The size of the magnetic fleld in the fluid rest frame in a flat 
space time is related to as 

\b\' = ^,B'+[v^B,]\ (2.60) 

Notice that there is a factor of An in difference with [18], due to different 
conventions for B^. 
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We may parametrise the problem, by the using the usual definition of 
P — ^/2P/\b\'^ as the ratio of gas to magnetic pressure in the fluid rest 
frame. For an ideal equation of state, in terms of /3 and P, ^ is written 




To facilitate comparison, I have used the same box size, < y < 3, pressure 
P = 1/3 10^^, background density p = I and amplitude of the perturbation, 
Aq = 10~^, as in [18]. The adiabatic index is relativistic with F = 4/3. The 
pulses are set up with a square formed wave in 

if 1 < 1/ < 1.5 

if 1.5 < y < 2 (2.62) 
elsewhere 

and for fixed p and P the Alfven velocities only depend on /? and w^. The 
parameters are given in table III and Fig. 2.11 shows at the time given in 
table III. The times have been selected to those moments in time where the 
two pulses line up exactly one after the other, and by visual inspection it is 

easy to see how the test fares. 

The amplitudes of the two waves are inversely proportional to their Lorentz 
factors [18] 

A- - W{vt) ^'-''^ 

and because the starting amplitude Aq is very small, the waves should not in- 
teract with each other. Then, the sum of the amplitudes is equal to the initial 
amplitude = A'^ + A~ . In table III, I have given the measured velocities 
and amplitudes together with the expected ones derived from Eqs. (2.59) and 
(2.63). 

The tests are selected to highlight different regimes of Eq. (2.59). In 
ALFl, we have a very low and consequently the Alfven velocity is close to 
the speed of light. The code does not fare well, showing 22% disagreement 
witht the expected value. In ALF2 the background fluid velocity is selected 
such that one pulse is frozen. It can be verifled from the flgure that the test 
is passed. In ALF3 v''^ — 0.8 and both pulses are travelling to the right. For 
the fast moving pulse, again the wave speed is too high with a 15% overshoot 
and the amplitudes are furthermore wrong. In ALF4 I have adjusted to 
yield two pulses with i>+ = "2^", and there are no problems with the test. 

The tests indicate that the code begins to signiflcantly overestimate the 
Alfen velocity when Va ^ 0.75, but in all cases the sum of the amplitudes is 
conserved. This is in accordance with the results from the shocktubes, where 
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correct jumps where observed albeit propagating with different velocities. 

The many tests presented in this section document both the strengths 
and weaknesses of the code. It is essential to know the limits of the code, 
not only in terms of stability, but also when to trust the physical models 
produced using it. 

It is clear that there are some stability problems with high Lorentz boosts, 
it is too viscous in the advection of high density blobs away from low density 
areas, and that it overestimates the Alfven speed, when it is relativistic. The 
cures to these problems are twofold: 

• The time derivatives of four velocities and the electric field have to be 
properly centred. 

• The time derivatives of the shear viscosity in Eqs. (2.21) and (2.22) 
have to be included. 

On the positive side the results all show flux conservation and reproduc- 
tion of the proper jump conditions across discontinuities both in HD and 
MHD tests. We can successfully model problems with severe pressure and 
density contrasts and in most cases faithfully resolve sharp features with 
very few points. This is done without showing oscillatory behaviour. Even 
though the largest fraction of the CPU time is spent calculating the shear 
viscosity, the gains in stability and the sharpness of discontinuous features 
increased fundamentally when I shifted from using a "mockup viscosity" to 
a full physically motivated one. 

The two points above are not fundamental or unsurmountable in any way 
and will be addressed in future work. 

2.5 Astrophysical Applications 

We can already apply the code to the understanding of mildly relativistic 
phenomena. Here I present first results from two applications related to the 
areas which motivated the development of the code. 

2.5.1 Decaying magnetic Gelds in the early universe 

In the introduction, we considered the evolution of magnetic fields in the 
early universe. Many analytical studies show that, at best, it will be very 
hard to find traces or fingerprints of primordial magnetic fields in the cosmic 
microwave background radiation, but these analyses do not take into account 
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Wave number k 



Fig. 2.12: Evolution of the power spectrum and the total energy density for a 



the non linear coupling between the different wave modes and the possibility 
of inverse cascades transferring energy from the small to the large scales. 

Christensson et al [15, 16] argued that, in fact, if a turbulent helical 
magnetic field was created, for example at the electro weak phase transition, 
it would undergo an inverse cascade, while a non- helical field would not. 

As a nontrivial 3D test of the code I have initialised a simple turbulent 
non-helical magnetic field and a turbulent velocity field with power spectra 
given as 



where the index k indicates the Fourier transform and due to causality, the 
exponents are constrained to > 0, > 2. In accordance with [15], I have 
taken them to be at their minimal value. The cut-off kc is introduced to limit 
numerical noise near the Nyquist frequency. In a 96^ run, with a box size of 
[0, 27r]^, where k = 1 corresponds to the largest mode in the box, I found that 
a value of k^ = 10 was sufficient to quench the numerical noise. To generate 
proper divergence free initial conditions, I first calculate the corresponding 
vector potential and then take the curl. The initial magnetic and kinetic 
energy are both 5 x 10~^ and the average density is p = 1. The internal 
energy is initialised such that the sound speed is relativistic, cj. = 1/3. 

Simulations of turbulence are very sensitive to the type and amount of 



turbulent magnetic field. The curves to the left are, in decaying order 
for t = [0,3,9,12,15,18]. 




(2.64) 




(2.65) 
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viscosity used, since it can alter the long time evolution of high frequency 
modes significantly. We have made a series of runs with less and less viscosity 
and noticed, that with coefficients ui^s < 0.003 and ^2 < 0.06 there was no 
change in the decay rate of the spectrum. To the left in Fig. 2.12 is shown 
the magnetic power spectrum at different times for a run with z/i 3 = 0.0003 
and z/2 = 0.006. Correspondingly, to the right is shown the evolution of the 
magnetic energy in the box that I find decays as Em oc Comparing my 
results with [15] I find good agreement. They found a Em oc t~^-^ scaling 
law. The main difference between my runs and theirs is that they evolve the 
non relativistic equations, while I use the relativistic equations. 

2.5.2 Relativistic jets 

Relativistic jets seem ubiquitous in the universe, to be found over a large 
range of scales, from sub parsecs to kilo parsecs [50]. It is one of the purest 
displays of special relativity at work. While the other application tested the 
impact of the viscosity in the code, a jet is an excellent test of the codes 
capability to handle strong shocks. 

Taking into account the large computer resources a 3D jet takes, I have 
chosen to make a 2D jet. At the moment, only Cartesian geometry is imple- 
mented in the code and I have constructed a slab jet, which is periodic in the 
2;-direction. The injection happens in the y-direction, where rigid boundaries 
are imposed, while there are periodic boundaries in the x-direction. To avoid 
significant collision of the bow shock with itself, the computational domain 
is a square with [N^, Ny) = (800, 800). I have tried both to inject the jet in 
an purely hydrodynamic medium void of magnetic fields and in a medium 
with a parallel magnetic field B^. As expected the main difference was fur- 
ther coUimation of the jet due to magnetic confinement. Similar experiments 
have been reported by other authors [41, 43]. The jet has an injection radius 
of Rj = 5 cells, the density contrast is PamUent/ pjet = 10 and the pressure 
is P = 1. There is pressure equilibrium between the jet and the ambient 
medium. We inject the jet with a Lorentz factor of 1^ = 1.5 and the adia- 
batic index is set to be relativistic with F = 4/3. In Fig. 2.13 a sequence of 
snapshots are shown. The large resolution and thin injection radius makes 
it possible to follow the jet until it becomes unstable and decays. At t = 500 
we see a classic jet. At the jet head there is a Mach shock, and material 
that has passed through the head is slowly forming a backfiow along the jet, 
building up a shear layer. Furthest out is the bow shock. The jet is unstable, 
and at later times the jet head disintegrates into a number of vortices and 
looses most of the kinetic energy. The perturbation runs backwards, slowly 
unwinding the spine of the jet. 




Fig. 2.13: From left to right: Density, pressure and four velocity. From top to 
bottom: The jet at t = 500, 1000, 1500, 2000. 
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2.6 Code Implementation: Performance and Scalability 

To successfully exploit modern massively parallel computers and clusters 
of computers, which at national centres of supercomputing often consist of 
hundreds of CPUs, the numerical implementation has to be carefully crafted. 
The program has to run at the optimal speed for small problems on a single 
CPU on a variety of architectures, while at the same time, it is important to 
distribute the workload evenly over all CPU's in the machine (be it a large 
shared memory computer or a cluster of off-the-shelf workstations). 



All of the fluid dynamics codes currently in use in Copenhagen are based on 
or derived from a common base code, the so called stagger code. The first ver- 
sion was made by Nordlund and Galsgaard in 1995 [61]. The GrMHD code 
makes use of the same basic principles. The equations of motion Eqs. (2.43), 
(2.44), (2.46), (2.47) are solved with finite differences through direct differen- 
tiation, and the variables are staggered on the grid. Scalar variables D and 
S, and derived scalar quantities are centred in each cell. The primary vector 
quantities Vi and are calculated on the faces of the cell while the derived 
vector quantities and J* are calculated on the edges (see Fig. 2.14). The 
boundary conditions are implemented as in [61]. 



Fig. 2.14: The basic staggering of different quantities on the grid. The figure was 
adapted from [24]. To make the figure visually easier to understand I 
have on purpose drawn a left handed coordinate system. I use a right 
handed coordinate system in the code. 

The differentiation operators are sixth order in space. Derivatives are 



2.6.1 The stagger code 
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calculated at half grid points and use a stencil of 6 points. In many cases 
this gives a natural placement of the variables, since the different quantities 
already are staggered in space. For example, the electric current J* is found 
located at the edges and according to Eq. (2.43) is the curl of (for the 
sake of simplicity in this example we disregard the displacement current, a 
and shear i3* is located at the face of each cell, and in the code the 
calculation of the current can be packed into three simple lines: 

r = ddydn(B^) - ddzdn(B^) 

= ddzdn(-B^) - ddxdn(S") (2.66) 
r = ddxdn(S^) - ddydn(6^) 

In some cases, most notably the complicated viscosity operator, the differen- 
tiation does not place the variables at the desired position on the grid, and 
interpolation has to be done. The corresponding interpolation operator is of 
fifth order. It also uses a stencil of 6 points [61]. A crucial addition to the 
original method described in [61], which later has been employed in most of 
the stagger based codes is the use of exponentials and logarithms to produce 
geometric interpolation. As an analogy, the geometric mean of two numbers 
may be rewritten in terms of the arithmetic mean of the logarithms 



G{a, b) = Va ■ b — exp 



- (In a + In b) 



= exp[ii'(lna,ln6)]. (2.67) 



Geometric interpolation has two very appealing qualities, when dealing with 
discontinuities across shocks. First of all, positive definite quantities, such as 
the density and the energy, stay positive. Secondly, geometric interpolation is 
a much better measure when the density or pressure is changing with orders 
of magnitude over a few points. This happens at shock fronts and surface 
transitions. 

To make the code easily readable and hide all the loops, where the dif- 
ferent interpolations and differentiations are done, all operators are hidden 
in a set of subroutines. In fact, Eq. (2.66) corresponds exactly to the simple 
version of the code. In the production version, we make an effort to optimise 
memory references and reuse the cache memory at each CPU, but even with 
full parallelisation the V x term only expands to 



Electric current I = curl(B) 



do kk=ks,ke 

call ddydn_set(Bz, Jx) ! Jx = ddydn(Bz) - ddzdii(By) 
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call ddzdn_sub(By , Jx) 

call ddxdn_set(By, Jz) ! Jz = ddxdn(By) - ddydn(Bx) 
call ddydn_sub(Bx, Jz) 

call ddzdn_set(Bx, Jy) ! Jy = ddzdii(Bx) - ddxdn(Bz) 
call ddxdii_sub(Bz, Jy) 
end do 

The reason why sixth order differentiation and fifth order interpolation 
operators are used in the stagger code is a question of balance between pre- 
cision and computational load. The highest wavenumber a given method 
can resolve depends not only on the resolution of the mesh, but also on the 
order of the interpolation and differentiation operators. It was found em- 
pirically by Nordlund and Galsgaard [61] that sixth order gives effectively a 
better resolution than fourth order operators, even after the somewhat larger 
computational cost of the 6th order operations are considered. Maron [48] 
made a formal investigation of the effective resolution of different methods 
and orders and found that a fourth order scheme can resolve up to 0.24 of 
the maximal wave number kmax, while a sixth order scheme resolves waves 
with up to 0.34 kmax- Going to eighth order the maximal wave number is 
0.4 kmax- At higher wave numbers the gain is negligible if one takes into 
account the added communication and amount of ghost zones that have to 
be allocated. 

2.6.2 The paper code: Optimal parallelisation and cache reuse 

Together with J. Hunsballe [37] I performed what essentially amounts to 
a complete rewrite of the stagger code. We still retain the qualities that 
have been described above. The basic physical equations are the same, and 
artificial viscosity is implemented in the same manner. We use the same 
high order interpolation and differentiation operators. The difference is in 
the technical details: Our goal has been to produce a very high level object 
oriented code, that is easily readable, runs at the highest possible speed and 
scale to hundreds of CPUs. 

In the stagger code the basic scope for any operator (such as a differen- 
tiation operator) has been the full three dimensional array. For example, to 
interpolate the density to face values one would write: 

call xdn_set(rho,xdnr) 
call ydn_set (rho ,ydnr) 
call zdn_set(rho,zdnr) 

The problem with this approach is that, on modern computers, the band- 
width between the main memory and the CPU is much lower than the com- 
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putational power. There is also a big latency involved. It can easily take 
200 clock cycles from the moment the CPU asks for a specific block of mem- 
ory until it actually is delivered. To alleviate this problem, there is a small 
amount of very fast memory, the cache, often integrated directly on the CPU. 
On current high end architectures, such as the Itanium, Power, Alpha, Sparc 
and Mips based machines the cache size is of the order of 3-8 MB per CPU, 
while a normal Opteron or Pentium based CPU has 1 MB cache. On a small 
problem, of todays standards, such as a 128^ mesh per CPU, the amount of 
cache taken by just one array is already 8 MB. That means the stagger code 
is already beginning to be memory bandwidth limited, and not limited by 
the speed of the CPU at this problem size. In the new Paper Code, the basic 
scope is instead a slice in the x — y plane (see Fig. 2.15); hence the name. 
The above hues of code would then be written: 

do kk=ks,ke 

call xdn_set (rhOjXdnr) 

call ydn_set (rho.ydnr) 

call zdn_set (rhOjZdnr) 
end 

where the kk-loop runs over the papers. The code is almost identical, since 
we hide the loop index in a global variable, but the characteristics are radi- 
cally different. Because we use a sixth order scheme we now only required to 
store 5 "papers" for the z-operator that needs values from different papers 
of p, in the cache to keep the CPU running at maximal speed. Even for 
a 512^ problem 5 papers only take up 5 MB of memory, and by testing on 
an SGI Altix machine with 3 MB of cache, we have found that performance 
starts to decrease around 400^, while at 1024^ x 20 performance has fallen 
to 2/3 of maximum. All modern CPUs are able to vectorise and pipeline 
simple instructions. By default, we have therefore chosen the innermost di- 
mension, the x-direction, to be as simple as possible, with periodic boundary 
conditions. Then, the compiler will be able to schedule essentially all opera- 
tions as SIMD instructions. The middle dimension, the y-direction, does not 
have any significance and is the best place to calculate boundary conditions, 
for problems that only contain boundaries in one direction. This way any 
computational load from the boundary is spread evenly over all the papers. 

So far in Copenhagen we have had good access to shared memory ma- 
chines. By far the easiest way to parallelise the code is then to use OpenMP. 
However, shared memory machines are relatively expensive and limited in 
size. A few versions exist of the stagger code that use MPI to run on clus- 
ters. One of the major current technology trends is the integration of two 
(Intel, AMD, Sun) or more (IBM) CPU cores on a single piece of sihcon. 
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BASIC COMPUTATIONAL BLOCK: A PAPER 




X Dimension: Innermost dimension in the loop. To keep it simple there are periodic 
boundaries and the CPU can use vector instructions such as SSE. 

Y Dimension: Middle dimension. This is the optimal place to apply boundary 
conditions. They are then distributed evenly among cpus. 

Z Dimension: Outer dimension. Used for parallelization. Basic operations are 
performed on one slice at constant z at a time. 



Fig. 2.15: The basic scope for any calculation in the new Paper Code is the x — y 
plane giving optimal reuse of cache, vectorisation and simple implemen- 
tation of boundary conditions. 



We can only expect that all CPUs in the future will be massively multi 
threaded. Then the optimal approach to parallelisation will be a hybrid one 
with OpenMP inside a single CPU node and MPI between nodes. Current 
and future parallelisation strategies have been sketched in Fig. 2.16. 

In the Paper Code we have effectively hidden the parallel nature of the 
code. Each CPU is automatically assigned a number of papers from ks to ke, 
and in the main part of the code, where the physics are calculated, one only 
has to consider dependencies in the z-direction and insert synchronisation 
points as appropriate. As an example consider the use of geometric means 
to interpolate the density to face values 

= exp(xdn(ln(p))) , py = exp(ydn(ln(p))) , = exp(zdn(ln(p))) , (2.68) 

where idn denotes interpolation half a point down in the i direction. This 
may be coded in two blocks: 

do kk=ks , ke 

lnr(:,:,kk) = alog(rho ( : , : ,kk) ) 
enddo 

!$omp barrier !< — Sync: zdnr=zdn_set (Inr) 

! depends on non-local papers 

do kk=ks , ke 

call xdn_exp(lnr ,xdnr) 
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PARALLELIZATION STRATEGIES 



CURRENT VERSION 












FUTURE VERSION 
Y: USING OPENMP 




Fig. 2.16: Parallelisation strategies. We have demonstrated perfect scalability up 
to 128 CPUs with our current OpenMP implementation. Future imple- 
mentations will be based on a hybrid OpenMP/MPI model. An added 
benefit of a hybrid model is the improved cache reuse for very large box 
sizes (i.e. 1024^ and beyond). 



call ydn_exp(lnr ,ydnr) 
call zdn_exp(lnr ,zdnr) 
end 

Notice that geometric interpolation is a common operation and to stream- 
line things we have made special interpolation operators, that automatically 
applies the exponential. 

A barrier works as a synchronisation point. The CPUs have to stop 
at a barrier and wait until all of them have arrived. When only using a 
small number of CPUs the number of barriers are not very important, but 
when considering hundreds of CPUs, it is essential that the barrier count 
is minimised. Any small disturbance for one CPU will make all the others 
wait at each barrier. To take an example: If there in each time step are 
100 barriers, and a CPU is randomly disturbed once during a time step, 
giving a slowdown of 1%, this extra noise will for two CPUS give rise to a 
2% slowdown. For hundred of CPUs the same disturbance, since it occurs 
in random sections, gives on average at least a 50% slowdown. By carefully 
analysing the numeric implementation, we have found that in a full update 
of the cells the minimum number of barriers needed to calculate any part of 
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Fig. 2.17: Scalability of the paper code: Results of scaling a simple HD experiment 
on an SGI Altix machine. To the left is the strong scalability for a 256^ 
experiment on a dedicated machine with a 1.3GHz CPU. To the right is 
shown the weak scalability running on a loaded machine with a 1.6GHz 
CPU with the experiment size varying according to the number of CPUs. 



the code is 6. The old stagger code based MHD was logically structured in 
different sections, according to the different physics. First, the calculation of 
velocities from momentum, then, the pressure, the viscosity, the stress tensor, 
the MHD equations and at last the equation of motion for the internal energy. 
Since all parts need between 4 and 6 barriers, one ends up having at least 
20 barriers. With the new code, we have applied a "principle of origami" 
folding the logical structure of the code. After each barrier, we consider 
all the different equations and calculate the maximum amount of physics. 
When threading the 5 small MHD parts, 6 viscosity parts etc together we 
end up having only 6 barriers. Recently, we had the opportunity to have 
the paper code tested on the NASA Columbia supercomputer and the code 
scaled efficiently on up to at least 128 CPUs. 

We have implemented a full HD/MHD code, including self gravity & 
turbulent driving, and a special relativistic HD/MHD version of the code 
described in this thesis. Both codes show spectacular performance. The 
MHD code can update 1.2 million cells per CPU per second and it runs at 
30% of theoretical peak performance on the SGI Altix machine. A normal 
grid based MHD code, even when optimised, performs at between 5% and 
10% of peak performance (see [62] for a detailed analysis of five state of the 
art codes). In fact, we believe that the paper code is one of the highest 
performing codes of its kind. This is both due to the low absolute cost of 
evaluating the MHD equations for a single cell and the effectiveness with 
which we have implemented the algorithm. The special relativistic MHD 
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version runs at 250.000 zone updates per second, which is also well above 
quoted numbers in the literature. 

2. 7 Discussion 

In this chapter I have discussed the theoretical foundation and numerical 
implementation of a new general relativistic magnetohydrodynamics code. 
When designing a new code without building upon existing work, it is tempt- 
ing to use an already existing theoretical basis. However, instead of this, I 
have derived a new form of the equations of motion with global coordinates 
evolving the dynamical variables from the point of view of a local observer. 
This approach makes it possible to employ a highly sophisticated artificial 
viscosity. This is just but an example of the possibilities the new formula- 
tion opens up for. The implication of my approach is that any new physics 
that is implemented and working in special relativity in the future, be it a 
new equation of state, radiative transfer, or a perturbative implementation 
of gravity, may easily be reused in the general relativistic version of the code. 
This may be done because the special and the general relativistic versions 
are related through the simple formulas given in App. B. When deriving the 
equations of motion, I have not made any assumptions about the background 
metric, so that the design is ready to be coupled with methods solving the 
full Einstein equations, such as the CactusCode. 

This new GrMHD code has been tested on a variety of demanding prob- 
lems, and it has been demonstrated that it is able to deal with huge pressure 
and density gradients. It shows some problems in the case of flows with high 
Lorentz factors, but they can be addressed and will be solved in the near 
future. The tests carried out include both synthetic benchmarks that tests a 
certain aspect of the code, and real astrophysical applications. 

The computer code is based on a refinement of the current infrastructure 
for fluid dynamics used in Copenhagen, which has been developed together 
with J. Hunsballe. It shows a spectacular performance on modern computer 
architectures, exploiting up to 30% of the theoretical peak performance. The 
special relativistic versions of the hydrodynamics and magnetohydrodynam- 
ics codes arc three dimensional and have been fully parallelised. They have 
been tested and scale to hundreds of CPUs, making it possible to exploit 
massive supercomputers at national centres to the full extent. 

I plan to employ the code in combination with the other numerical tools 
presented in this thesis in order to understand extreme astrophysics near 
compact objects. A first joint application of the particle code and this code 
is presented in Chapter 5. Furthermore, observational cosmology is reaching 
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a level of quality, where soon not everything can be addressed in terms of 
simple one dimensional linear perturbation theory, and I plan to employ the 
code in the understanding of the non-trivial evolution of magnetic fields in 
the early universe. 
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3. MAGNETIC FIELD GENERATION IN COLLISIONLESS 
SHOCKS; PATTERN GROWTH AND TRANSPORT 



In this chapter I present results from three-dimensional particle simulations of 
coUisionless shock formation, with relativistic counter-streaming ion-electron 
plasmas first published in Fredriksen et al. [26]. Particles are followed over 
many skin depths downstream of the shock. Open boundaries allow the 
experiments to be continued for several particle crossing times. The exper- 
iments confirm the generation of strong magnetic and electric fields by a 
Weibel-like kinetic streaming instability, and demonstrate that the electro- 
magnetic fields propagate far downstream of the shock. The magnetic fields 
are predominantly transversal, and are associated with merging ion current 
channels. The total magnetic energy grows as the ion channels merge, and 
as the magnetic field patterns propagate down stream. The electron pop- 
ulations are quickly thermalised, while the ion populations retain distinct 
bulk speeds in shielded ion channels and thermalise much more slowly. The 
results help us to reveal processes of importance in coUisionless shocks, and 
may help to explain the origin of the magnetic fields responsible for afterglow 
synchrotron/jitter radiation from Gamma- Ray Bursts. 

3.1 Introduction 

The existence of a strong magnetic field in the shocked external medium is 
required in order to explain the observed radiation in Gamma- Ray Burst 
afterglows as synchrotron radiation [e.g. 65]. Nearly coUisionless shocks, 
with synchrotron-type radiation present, are also common in many other 
astrophysical contexts, such as in supernova shocks, and in jets from active 
galactic nuclei. At least in the context of Gamma-Ray Burst afterglows the 
observed synchrotron radiation requires the presence of a stronger magnetic 
field than can easily be explained by just compression of a magnetic field 
already present in the external medium. 

Medvcdcv & Loeb [53] showed through a linear kinetic treatment how a 
two-stream magnetic instability - a generalisation of the Weibel instability 
[81, 84] - can generate a strong magnetic field (eg, defined as the ratio of 
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magnetic energy to total kinetic energy, is 10~^-10~^ of equipartition value) 
in coUisionless shock fronts [see also discussion in 71]. We note in passing 
that this instability is well-known in other plasma physics disciplines, e.g. 
laser-plasma interactions [14, 82], and has been applied in the context of 
pulsar winds by Kazimura et al. [38]. 

Using three-dimensional particle-in-cell simulations to study relativistic 
coUisionless shocks (where an external plasma impacts the shock region with 
a bulk Lorentz factor F = 5 — 10), Frederiksen et al. [25], Nishikawa et al. [60], 
and Silva et al. [73] investigated the generation of magnetic fields by the two- 
stream instability. In these first studies the growth of the transverse scales 
of the magnetic field was limited by the dimensions of the computational 
domains. The durations of the Nishikawa et al. [60] experiments were less 
than particle travel times through the experiments, while Silva et al. [73] 
used periodic boundary conditions in the direction of streaming. Further, 
Frederiksen et al. [25] and Nishikawa et al. [60] used electron-ion {e~p) plas- 
mas, while experiments reported upon by Silva et al. [73] were done with 
e~e"'" pair plasmas. 

Here, we report on 3D particle-in-cell simulations of relativistically counter- 
streaming e'p plasmas. Open boundaries are used in the streaming direction, 
and experiment durations are several particle crossing times. Our results can 
help to reveal the most important processes in coUisionless shocks, and help 
to explain the observed afterglow synchrotron radiation from Gamma-Ray 
Bursts. We focus on the earhest development in shock formation and field 
generation. Late stages in shock formation will be addressed in Chapter 5. 

3.2 Simulations 

Experiments were performed using a self-consistent 3D3V electromagnetic 
particle-in-cell code originally developed for simulating reconnection topolo- 
gies [36] , redeveloped by Frederiksen [24] to obey special relativity and to be 
second order accurate in both space and time. 

The code solves Maxwell's equations for the electromagnetic field with 
continuous sources, with fields and field source terms defined on a stag- 
gered 3D Yee-latticc [83]. The sources in Maxwell's equations are formed 
by weighted averaging of particle data to the field grid, using quadratic 
spline interpolation. Particle velocities and positions are defined in con- 
tinuous (r, 7v)-space, and particles obey the relativistic equations of motion. 

The grid size used in the main experiment was (x, y, z) = 200 x 200 x 
800, with 25 particles per cell, for a total of 8 x 10^ particles, with ion to 
electron mass ratio mi /me — 16. To adequately resolve a significant number 
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Fig. 3.1: The left hand side panel shows the longitudinal electron current density 
through a transverse cut at z = 100, with a small inset showing the 
ion current in the same plane. The right hand side panel shows the ion 
current at z = 600 = 30(5j, with the small inset now instead showing 
the electron current. The arrows represent the transverse magnetic field. 
Both panels are from time t = 1200. 



of electron and ion skin-depths {Se and 6i), the box size was chosen such that 
L^^y = I05i ~ iOSe and ~ 405j ~ 1606 e. Varying aspect and mass ratios 
were used in complementary experiments. 

Two counter-streaming - initially quasi-neutral and cold - plasma popu- 
lations are simulated. At the two-stream interface (smoothed around z = 80) 
a plasma {z < 80) streaming in the positive z-direction, with a bulk Lorentz 
factor r = 3, hits another plasma {z > 80) at rest in our reference frame. 
The latter plasma is denser than the former by a factor of 3. Experiments 
have been run with both initially sharp and initially smooth transitions, 
with essentially the same results. The long simulation time gradually al- 
lows the shock to converge towards self-consistent jump conditions. Periodic 
boundaries are imposed in the x- and y-directions, while the boundaries at 
z = and z = 800 are open, with layers absorbing transverse electromagnetic 
waves. Inflow conditions at 2; = are fixed, with incoming particles supplied 
at a constant rate and with uniform speed. At z = 800 there is free outflow 
of particles. The maximum experiment duration is 480 uj~^ (where is 
the electron plasma frequency), sufficient for propagating F ^ 3 particles 2.8 
times through the box. 
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3.3 Results and Discussions 

The extended size and duration of these experiments make it possible to 
follow the two-stream instability through several stages of development; first 
exponential growth, then non-linear saturation, followed by pattern growth 
and downstream advection. We identify the mechanisms responsible for these 
stages below. 

3.3.1 Magnetic Geld generation, pattern growth 
and Held transport 

Encountering the shock front the incoming electrons are rapidly (being lighter 
than the ions) deflected by flcld fluctuations growing due to the two-stream 
instability [53]. The initial perturbations grow non-linear as the deflected 
electrons collect into first caustic surfaces and then current channels (Fig. 3.1). 
Both streaming and rest frame electrons are defiected, by arguments of sym- 
metry. 

In accordance with Ampere's law the current channels are surrounded by 
approximately cylindrical magnetic fields (illustrated by arrows in Fig. 3.1), 
causing mutual attraction between the current channels. The current chan- 
nels thus merge in a race where larger electron channels consume smaller, 
neighbouring channels. In this manner, the transverse magnetic field grows 
in strength and scale downstream. This continues until the fields grow strong 
enough to defiect the much heavier ions into the magnetic voids between the 
electron channels. The ion channels are then subjected to the same growth 
mechanism as the electrons. When ion channels grow sufficiently powerful, 
they begin to experience Debye shielding by the electrons, which by then 
have been significantly heated by scattering on the increasing electromag- 
netic field structures. The two electron populations, initially separated in 
7v-space, merge to a single population in approximately 205e {z = 80-200) 
as seen in Fig. 3.6. The same trend is seen for the ions - albeit the merg- 
ing rate might be significantly slower than predicted by extrapolating with 
mi/rrie, since Debye shielding stabilises the ion channels. 

The Debye shielding quenches the electron channels, while at the same 
time supporting the ion-channels; the large random velocities of the electron 
population allow the concentrated ion channels to keep sustaining strong 
magnetic fields. Fig. 3.1, shows the highly concentrated ion currents, the 
more diffuse - and shielding - electron currents, and the resulting magnetic 
field. The electron and ion channels are further illustrated in Fig. 3.2. Note 
the limited 2;-extent of the electron current channels, while the ion current 
channels extend throughout the length of the box, merging to form larger 
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Fig. 3.2: Electron (top) and ion (bottom) currents, averaged over the x-direction, 
at time t = 1200. 



scales downstream. Because of the longitudinal current channels the mag- 
netic field is predominantly transversal; we find ~ 10~^ — 10~^. 

Figure 3.3 shows the temporal development of the transverse magnetic 
field scales around z = 250. The power spectra follow power- laws, with the 
largest scales growing with time. The dominant scales at these z are of the 
order 6i at early times. Later they become comparable to L^^y. Figure 3.4 
captures this scaling behaviour as a function of depth for t = 2400. 

The time evolutions of the electric and magnetic field energies are shown 
in Fig. 3.5. Seeded by fluctuations in the fields, mass and charge density, the 
two-stream instability initially grows super-linearly (t = 80 — 100), refiecting 
approximate exponential growth in a small sub-volume. Subsequently the 
total magnetic energy grows more linearly, refiecting essentially the increasing 
volume filling factor as the non-linearly saturated magnetic field structures 
are advected downstream. 

At t ~ 1100 the slope drops off, due to advection of the generated fields 
out of the box. The continued slow growth, for t > 1100, reflects the increase 
of the pattern size with time (cf. Fig. 3.3). A larger pattern size corresponds, 
on the average, to a larger mean magnetic energy, since the total electric 
current is split up into fewer but stronger ion current channels. The magnetic 
energy scales with the square of the electric current, which in turn grows in 
inverse proportion to the number of current channels. The net effect is that 
the mean magnetic energy increases accordingly. 
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Fig. 3.3: Power spectrum of for z = 250 at different times. 



The magnetic energy density keeps growing throughout our experiment, 
even though the duration of the experiment (480 cOpJ-) significantly exceeds 
the particle crossing time, and also exceeds the advection time of the mag- 
netic field structures through the box. This is in contrast to the results 
reported by Silva et al. [73], where the magnetic energy density drops back 
after about 10-30 Up^. It is indeed obvious from the preceding discussion 
that the ion-electron asymmetry is essential for the survival of the current 
channels. 

From the requirement that the total plasma momentum should be con- 
served, the (electro) magnetic field produced by the two-stream instability 
acquires part of the z-momentum lost by the two-stream population in the 
shock; this introduces the possibility that magnetic field structures created 
in the shock migrate downstream of the shock and thus carry away some of 
the momentum impinging on the shock. 

Our experiments show that this does indeed happen; the continuous in- 
jection of momentum transports the generated field structures downstream 
at an accelerated advection speed. The dragging of field structures through 
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Fig. 3.4: Relative electromagnetic energy density cb. The contour colour plot 
shows the power in the transverse magnetic field through the box dis- 
tributed on spatial Fourier modes at t = 2400, with the dotted line mark- 
ing the wavenumber with maximum power. Superposed is the spatial dis- 
tribution of e^, averaged across the beam, at t = 2320 (dashed-dotted) 
and t = 2400 (full drawn), highlighting how EM- fields are advected down 
through the box. 



the dense plasma acts as to transfer momentum between the in-streaming 
and the shocked plasmas. 

3.3.2 Thermalisation and plasma heating 

At late times the entering electrons are effectively scattered and thermalised: 
The magnetic field isotropises the velocity distribution whereas the electric 
field generated by the e~-p charge separation acts to thermalise the popula- 
tions. Figure 3.6 shows that this happens over the ~ 20 electron skin depths 
from around z = 80 - 200. The ions are expected to also thermalise, given 
sufficient space and time. This fact leaves the massive ion bulk momentum 
constituting a vast energy reservoir for further electron heating and acceler- 
ation. Also seen in Fig. 3.6, the ions beams stay clearly separated in phase 
space, and are only slowly broadened (and heated). 

We do not see indications of a super-thermal tail in the heated electron 
distributions, and there is thus no sign of second order Fermi-acceleration 
in the experiment presented in this Letter. [60] and [73] reported accelera- 
tion of particles in experiments similar to the current experiment, except for 
more limited sizes and durations, and the use of an e~e~^ plasma [73]. On 
closer examination of the published results it appears that there is no actual 
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Fig. 3.5: Total magnetic (full drawn) and electric (dashed) energy in the box as a 
function of time. The inset shows a log-log plot of the same data. 



disagreement regarding the absence of accelerated particles. Whence, [60] 
refer to transversal velocities of the order of 0.2c (their Fig. 3b), at a time 
where our experiment shows similar transversal velocities (cf. Fig. 3.6) that 
later develop a purely thermal spectrum. [73] refer to transversal velocity 
amplitudes up to about 0.8c (their Fig. 4), or vy ~ 2, with a shape of the 
distribution function that appears to be compatible with thermal. In com- 
parison, the electron distribution illustrated by the scatter plot in Fig. 3.6 
covers a similar interval of v^, with distribution functions that are close 
to Lorentz boosted relativistic MaxweUians (see App. A for a discussion of 
Lorentz boosted thermal profiles). Thus, in the experiment reported on in 
this chapter there is no compelling evidence for non-thermal particle accel- 
eration. Thermalisation is a more likely cause of the increases in transversal 
velocities. 

Frederiksen et al. [25] reported evidence for particle acceleration, with 
electron gammas up to ~ 100, in experiments with an external magnetic field 
present in the up-stream plasma. This is indeed a more promising scenario for 
particle acceleration experiments (although in the experiments by [60] results 
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Fig. 3.6: Thermalisation and longitudinal acceleration, illustrated by scatter plots 
of the electron (orange) and ion (blue) populations. Note the back- 
scattered electron population {vz'jiv) < 0). 



with an external magnetic field were similar to those without). Figure 3.6 
shows the presence of a population of back-scattered electrons {vz'y < 0). In 
the presence of an external magnetic field in the in-streaming plasma, this 
possibly facilitates Fermi acceleration in the shock. 

3.4 Conclusions 

The experiment reported upon in this chapter illustrates a number of funda- 
mental properties of relativistic, collisionless shocks: 

1. Even in the absence of a magnetic field in the up-stream plasma, 
a small scale, fluctuating, and predominantly transversal magnetic field is 
unavoidably generated by a two-stream instability reminiscent of the Weibel- 
instability. In the current experiment the magnetic energy density reaches a 
few percent of the energy density of the in-coming beam. 

2. In the case of an e~p plasma the electrons are rapidly thermalised. 
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while the ions form current channels that are the sources of deeply pene- 
trating magnetic field structures. The channels merge in the downstream 
direction, with a corresponding increase of the average magnetic energy with 
shock depth. This is expected to continue as long as a surplus of bulk relative 
momentum remains in the counter-streaming plasmas. 

3. The generated magnetic field patterns are advected downstream at 
speeds intermediate of the streaming and rest frame plasmas. The electro- 
magnetic field structures thus provide scattering centres that interact with 
both the fast, in-coming plasma, and with the plasma that is initially at rest. 
As a result the electron populations of both components quickly thermalise 
and form a single, Lorentz-boosted thermal electron population. The two 
ion populations merge much more slowly, with only gradually increasing ion 
temperatures. 

4. The observed strong turbulence in the field structures at the shocked 
streaming interface provides a promising environment for particle accelera- 
tion. 

We emphasise that quantification of the interdependence and develop- 
ment of eu and es is accessible by means of such experiments as reported 
upon here. 

Rather than devising abstract scalar parameters and eu, that may be 
expected to depend on shock depth, media densities etc., a better approach 
is to compute synthetic radiation spectra directly from the models, and then 
apply scaling laws to predict what would be observed from corresponding, 
real supernova remnants and Gamma-Ray Burst afterglow shocks. 



4. NON-FERMI POWER LAW ACCELERATION IN 
ASTROPHYSICAL PLASMA SHOCKS 



CoUisionless plasma shock theory, which apphes for example to the afterglow 
of gamma ray bursts, still contains key issues that are poorly understood. In 
this chapter I discuss the results of charged particle dynamics in a highly rela- 
tivistic coUisionless shock numerically using ~ 10^ particles first published by 
Hededal et al. [32] . We find a power law distribution of accelerated electrons, 
which upon detailed investigation turns out to originate from an acceleration 
mechanism that is decidedly different from Fermi acceleration. Electrons are 
accelerated by strong filamentation instabilities in the shocked interpene- 
trating plasmas and coincide spatially with the powerlaw distributed current 
filamentary structures. These structures are an inevitable consequence of 
the now well established Weibel-like two-stream instability that operates in 
relativistic coUisionless shocks. The electrons are accelerated and deceler- 
ated instantaneously and locally; a scenery that differs qualitatively from 
recursive acceleration mechanisms such as Fermi acceleration. The slopes 
of the electron distribution powerlaws are in concordance with the particle 
powerlaw spectra inferred from observed afterglow synchrotron radiation in 
gamma ray bursts, and the mechanism can possibly explain more generally 
the origin of non-thermal radiation from shocked inter- and circum-stellar 
regions and from relativistic jets. 



4. 1 Introduction 



Given the highly relativistic conditions in the outflow from gamma ray bursts 
(GRBs), the mean free path for particle Coulomb collisions in the afterglow 
shock is several orders of magnitude larger than the fireball itself. In ex- 
plaining the microphysical processes that work to define the shock, MHD 
becomes inadequate and coUisionless plasma shock theory stands impera- 
tive. In particular two key issues remain, namely the origin and nature of 
the magnetic field in the shocked region, and the mechanism by which elec- 
trons are accelerated from a thermal population to a powerlaw distribution 
N{'y)d^ (X 7~^. Both ingredients are needed to explain the observed after- 
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glow spectra (e.g. [45, 64]). 

Regarding the origin of the magnetic field in the shocked region, obser- 
vations are not compatible with a compressed inter-stellar magnetic field, 
which would be orders of magnitude smaller than needed [30]. It has been 
suggested that a Weibel-like two-stream instability can generate a magnetic 
field in the shocked region (see Chapter 3, and Medvedev & Loeb [53]; Fred- 
eriksen et al. [25]; Nishikawa et al. [60]; Silva et al. [73]). Computer experi- 
ments presented in Chapter 3 and [26] showed that the nonlinear stage of a 
two-stream instability induces a magnetic field in situ with an energy con- 
tent of a few percent of the equipartition value, consistent with that required 
by observations. 

Fermi acceleration [22] has, so far, been widely accepted as the mechanism 
that provides the inferred electron acceleration. It has been employed exten- 
sively in Monte Carlo simulations (e.g. [59] and references therein), where it 
operates in conjunction with certain assumptions about the scattering of par- 
ticles and the structure of the magnetic field. The mechanism has, however, 
not been conclusively demonstrated to occur in ab initio particle simulations. 
As pointed out by Niemiec & Ostrowski [59], further significant advance in 
the study of relativistic shock particle acceleration is unlikely without under- 
standing the detailed microphysics of collisionlcss shocks. Also, recently Bar- 
ing and Braby [7] found that particle distribution functions (PDFs) inferred 
from GRB observations are in contradistinction with standard acceleration 
mechanisms such as diffusive Fermi acceleration. 

In this chapter we study ab initio the particle dynamics in a coUisionless 
shock with bulk Lorentz factor F = 15. We find a new particle acceleration 
mechanism, which is presented in section 4.2. Detailed numerical results 
are presented and interpreted in section 4.3, while section 4.4 contains the 
conclusions. 

4.2 A New Acceleration Mechanism 

A series of numerical experiments have been performed where coUisionless 
shocks are created by two colliding plasma populations. These experiments 
are described in more detail below, but a common feature is that the elec- 
tron PDF has a high energy tail which is powerlaw distributed. By carefully 
examining the paths of representative accelerated electrons, tracing them 
backwards and forwards in time, it has been possible to identify the mecha- 
nism responsible for their acceleration. The acceleration mechanism, which 
was presented for the first time in [32], works as follows: 

When two non-magnetised coUisionless plasma populations interpene- 
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Fig. 4.1: (A) Ray traced electron paths (red) and current density (blue). The 
colours of the electron paths reflect their four velocity according to the 
colour table in the inset (B). The shadows are equivalent to the x and 
y projections of their paths. The ion current density is shown with blue 
colours according to the colour table in the inset. The inset also shows 
the ion current density (blue) integrated along the x axis with the spatial 
distribution of fast moving electrons (red) over plotted. 



trate, current channels are formed through a Weibel-hke two-stream instabil- 
ity (see Chapter 3; Medvedev & Loeb [53]; Frederiksen et al. [25]; Nishikawa 
et al. [60]; Silva et al. [73]). In the nonlinear stage of evolution of this insta- 
bility, ion current channels merge into increasingly stronger patterns, while 
electrons act to Debye shield these channels, as shown in Chapter 3. Further 
it was showed that a Fourier decomposition of the transverse structure of the 
ion current filaments exhibits powerlaw behaviour which has been recently 
confirmed by Medvedev et al. [52]. 

At distances less than the Debye length, the ion current channels are sur- 
rounded by transverse electric fields that accelerate the electrons toward the 
current channels. However, the magnetic fields that are induced around the 
current channels act to deflect the path of the accelerated electrons, boost- 
ing them instead in the direction of the ion flow. Since the forces working 
are due to quasi-stationary fields the acceleration is a simple consequence of 
potential energy being converted into kinetic energy. Therefore the electrons 
are decelerated again when leaving the current channel, and reach their max- 
imal velocities at the centres of the current channels. Hence, as illustrated 
by Fig. 4. IB, the spatial distribution of the high energy electrons is a di- 
rect match to the ion current channels and the properties of the accelerated 
electrons depend primarily on the local conditions in the plasma. 
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Fig. 4.2: An ion current channel surrounded by an electric - and a magnetic field. 

Electrons in the vicinity of the current channels are thus subject to a 
Lorentz force with both an electric and magnetic component, working 
together to accelerate the electrons along the ion flow. Crossing the centre 
of the channel the process reverses leading to an oscillating movement 
along the channel. 



One might argue that the near-potential behaviour of the electrons, where 
they essentially must loose most of their energy to escape from the current 
channels, would make the mechanism uninteresting as an acceleration mech- 
anism since fast electrons cannot easily escape. However, this feature may 
instead be a major advantage, since it means that energy losses due to escape 
are small, and that the electrons remain trapped long enough to have time to 
loose their energy via a combination of bremsstrahlung and synchrotron or 
jitter radiation. We observe that only a very small fraction of the electrons 
manage to escape, while still retaining most of their kinetic energy. This 
happens mainly at sudden bends or mergers of the ion channels, where the 
electron orbits cannot be described in terms of a particle moving in a static 
electromagnetic field. 

To analyse the acceleration scenario quantitatively we construct a toy 
model. It has been sketched in Fig. 4.2. We assume that the ion current 
channel has radius R, that the total charge inside the cylinder per unit length 
is A and the ions all stream with velocity u and Lorentz factor F in the 
laboratory rest frame (see Fig. 4.2 and inset for definition of rest frames). 

Consider an electron with charge —q and mass m at a distance r from the 
centre of the channel, initially having no velocity components perpendicular 
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to the cylinder, and four velocity 70^^,0 parallel to the cylinder, and disregard 
for the moment any other shielding electrons. 

By analysing everything in the ion channel rest frame the problem reduces 
to electrostatics and it is possible to analytically calculate the change in four 
velocity of the electron when it reaches the surface of the cylinder. In the 
ion channel rest frame the electron has the Lorentz factor and four velocity 

^',^r^o(l-uv,,o), (4.1) 
7o^U = r7o(v,,o - u) , (4.2) 

where quantities in the ion channel rest frame are denoted with a prime. 
The ions were before moving with velocity u, and hence subject to a Lorentz 
contraction, but are now in their rest frame. The line charge density is 
reduced by a factor of F: A' = A/F. The electron will be attracted to 
the cylinder and will gain downward momentum in the r-direction. This is 
simply a free fall in the electric potential and the final velocity, when the 
electron reaches the outer edge of the cylinder, can be found by calculating 
the change in potential energy 



^K^t"" = - / dr = Hr/R) . (4.3) 



R 27reo 

The change in the Lorentz boost 7' is then mc^ A7' = ^E'f^^^ = —AEp^^. 
The electric force only works along the r-axis and the four velocity along the 
z-Sixis of the electron is conserved in the ion channel rest frame. Exploiting 
this we can calculate not only the total change in energy but also the change 
in the different velocity components. Returning to the laboratory rest frame 
we find 

^electron = 7^ 5" — , (4.4) 

^ilVz) electron = uA^electron ■ (4.5) 

The change in the Lorentz boost is directly proportional to the total charge 
inside the channel and inversely proportional to the electron mass. In reality 
the Debye shielding reduces the electric field further away from the ion chan- 
nel, so the estimate above is only valid for distances smaller than a Debye 
length. Inside the ion channel the electron is accelerated too, but the amount 
depends on the detailed charge distribution of the ions, and one should re- 
member, that in general the electrons do indeed have velocity components 
perpendicular. The above estimate then can be understood as an upper limit 
to the observed acceleration. 
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4.3 Computer Experiments 

The experiments were performed with the three-dimensional relativistic ki- 
netic and electromagnetic particle-in-ccll code described briefly in Chapter 3 
and more thoroughly in [24]. The code works from first principles, by solving 
Maxwell's equations for the electromagnetic fields and solving the Lorentz 
force equation of motion for the particles. 

Two colliding plasma populations are set up in the rest frame of one of 
the populations (downstream, e.g. a jet). A less dense population (upstream, 
e.g. the ISM) is continuously injected at the left boundary with a relativistic 
velocity corresponding to a Lorentz factor F = 15. The two populations 
initially differ in density by a factor of 3. We use a computational box with 
125 X 125 X 2000 grid points and a total of 8 x 10^ particles. The ion rest 
frame plasma frequency in the downstream medium is oOpi = 0.075, rendering 
the box 150 ion skin depths long. The electron rest frame plasma frequency 
is ujpe = 0.3 in order to resolve also the microphysics of the electrons. Hence 
the ion-to-electron mass ratio is mi/rrie — 16. Other mass ratios and plasma 
frequencies were used in complementary experiments. Initially, both plasma 
populations are unmagnetised. 

The maximum experiment duration has tmax = 340 cj~/, which is suf- 
ficient for the continuously injected upstream plasma (F = 15, ~ c) to 
travel 2.3 times the length of the box. The extended size and duration 
of these experiments enable observations of the streaming instabilities and 
concurrent particle acceleration through several stages of development [26]. 
Momentum losses to radiation (cooling) are presently not included in the 
model. We have, however, verified that none of the accelerated particles 
in the experiment would be subject to significant synchrotron cooling. The 
emitted radiation may thus be expected to accurately reflect the distribution 
of accelerated electrons. 

When comparing numerical data with Eq. (4.4) we take r to be the ra- 
dius where Debye shielding starts to be important. Using a cross section 
approximately in the middle of Fig. 4.1 we find A^jVz) electron = 581n(r/i?). 
It is hard to determine exactly when Debye shielding becomes effective, but 
looking at electron paths and the profile of the electric field we estimate that 
ln{r/R) ^ 1.3. Consequently, according to Eq. (4.4), the maximally attain- 
able four velocity in this experiment is in the neighbourhood of {'jVz)max = 75. 
This is in good agreement with the results from our experiments, where the 
maximum four velocity is {'yVz)max — 80. 

The theoretical model does of course not cover all details of the experi- 
ment. For example, in general the electrons also have velocity components 
parallel to the magnetic field; instead of making one dimensional harmonic 
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Fig. 4.3: A scatter plot of the local ion current density Jjon versus the four velocity 
of the electrons in a region downstream of the shock. Over plotted is a 
line (thin) showing the average four velocity as a function of J/on; and 
a line (thick) showing a straight line fit. Because 'cold' trapped thermal 
electrons (indicated with the ellipse) exist inside the ion current channel 
they count towards lowering the average four velocity at high Jjon- If 
the scatter plot was cleaned, statistically removing all thermal electrons, 
we would see a much tighter relation. Such cleaning, though, is rather 
delicate and could introduce biases by itself. The trend is clearly there 
though even for the 'raw' data. 



oscillations in the plane perpendicular to the current channel the electrons 
will describe complicated ellipsoidal paths. Fig. 4.1A shows the path of two 
electrons in the vicinity of an ion channel. But, overall, the electrons be- 
have as expected from the model considerations. Consequently, high speed 
electrons are tightly coupled to the ion channels, as clearly illustrated by 
Fig. 4.1B. 

Figure 4.4 shows that the electrons are powerlaw distributed at high ener- 
gies, with index p = 2.7. The electrons at the high gamma cut-off are found 
where the ion current peaks, as may be seen from Fig. 4.3. The maximum 
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ion current is limited by the size of our box; larger values would probably 
be found if the merging of current channels could be followed further down 
stream. The PDF is not isotropic in any frame of reference due to the high 
anisotropy of the Weibel generated electromagnetic field. The powerlaw in 
the electron PDF is dominant for 10 < 7 < 30. Likewise, a powerlaw dom- 
inates the ion current channel strength, Jjon, for 100 < J/on < 1000 (inset). 
A relation between the powerlaw distributions of these two quantities on 
their respective intervals is provided with Fig.4.3: We see that the aver- 
age four velocity is proportional (straight fine fit) to a power of the local 
ion current density on their respective relevant intervals, 10 < 7 < 30 and 
100 < Jion < 1000. Their kinship stems from the fact that acceleration is 
local. J/on has a powerlaw tail and its potential drives the high energy dis- 
tribution of the electrons according to Eq. (4.4), thus forming a powerlaw 
distributed electron PDF. 

Measuring the rate at which the in-streaming ions transfer momentum to 
the ion population initially at rest allows us to make a crude estimate of the 
length scales over which the two-stream instability in the current experiment 
would saturate due to ion thermalisation. A reasonable estimate appears to 
be approximately 10 times the length of the current computational box, or 
about 1500 ion skin depths. Assuming that the shock propagates in an in- 
terstellar environment with a plasma density of ~ 10® m~^ we may calculate 
a typical ion skin depth. Comparing this value with the upstream ion skin 
depth from our experiments, we find that the computational box corresponds 
to a scale of the order of 10^ m, or equivalently that the coUisionless shock 
transition region of the current experiment corresponds to about 10^ m. For 
an ion with a Lorentz factor 7 = 15 this length corresponds roughly to 40 
ion gyro radii in the average strength of the generated magnetic field. But it 
should be stressed that the in-streaming ions actually do not really gyrate 
since they mainly travel inside the ion current channels where the magnetic 
field, by symmetry, is close to zero. Also, the strong electromagnetic fields 
generated by the Weibel instability and the non-thermal electron accelera- 
tion, which is crucial from the interpretation of GRB afterglow observations, 
emphasise the shortcoming of MHD in the context of coUisionless shocks. 

In the computer experiments presented here we have used a mass ratio 
mi/rrif. = 16 in order to resolve the dynamics of both species. Eq. (4.4) 
suggests that reducing the electron mass to 1/1836 will increase the ac- 
celeration of the electrons, but the gained energy is independent of the mass 
(see Eq. (4.3)). In this experiment we observe electrons with energies of 
approximately 5 GeV. Even further acceleration may occur as ion channels 
keep growing down stream, outside of our computational box. 

The scaling estimates above depend, among other things, on plasma den- 
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sities, the bulk Lorcntz factor, and the mass ratio (mi/nie). A parameter 
study is necessary to explore these dependencies, but this is beyond the 
scope of the present chapter. We thus stress that the extrapolations per- 
formed here are speculative and that unresolved physics could influence the 
late stages of the instability in new and interesting ways as discussed in the 
following chapter. 

When the in-streaming ions are fully thermalised they can no longer 
support the magnetic field structures. Thus one might speculate that the 
radiating region of the GRB afterglow is actually very thin, as suggested by 
Rossi & Rees [71]. Further, traditional synchrotron radiation theory does not 
apply to intermittent magnetic field generated by the two-stream instability, 
since the electron gyro radii often are larger than the scales of the magnetic 
field structures. We emphasise the importance of the theory of jitter radiation 
for understanding the generated radiation [51]. 

4.4 Conclusions 

In this chapter we have proposed an acceleration mechanism for electrons 
in coUisionless shocks. The theoretical considerations were suggested by 
particle-in-cell computer experiments, which also allowed quantitative com- 
parisons with the theoretical predictions. We have shown that the non- 
thermal acceleration of electrons is directly related to the ion current chan- 
nels in the shock transition zone. The results are applicable to interactions 
between relativistic outflows and the interstellar medium. Such relativistic 
outflows occur in GRB afterglows and in jets from compact objects [21]. 
The suggested acceleration scenario might overcome some of the problems 
pointed out by Baring & Braby [7] regarding the apparent contradiction be- 
tween standard Fermi acceleration and spectral observations of GRBs. 

The mechanism has important implications for the way we understand 
and interpret observations of coUisionless shocks: 

1 . The acceleration mechanism is capable of creating a powerlaw electron 
distribution in a coUisionless shocked region. In the computer experiment 
presented here a bulk flow with F = 15 results in a powerlaw slope p = 2.7 
for the electron PDF. Additional experiments will be needed to disentangle 
what determines the exact value of the slope. 

2. The acceleration is local; electrons arc accelerated to a powerlaw in 
situ. Therefore the observed radiation field may be tied directly to the lo- 
cal conditions of the plasma and could be a strong handle on the physical 
processes. 

3. Our results strengthen the point already made in Chapter 3; that the 
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Fig. 4.4: The normalised electron particle distribution function downstream of the 
shock. The dot-dashed line is a powerlaw fit to the non-thermal high 
energy tail, while the dashed curve is a Lorentz boosted thermal electron 
population. The histogram is made from the four velocities of electrons 
in a thin slice in the z-direction of the computational box. The inset 
shows a similar histogram for ion current density sampled in each grid 
point in the same slice. The bump in the inset is a statistical fluctuation 
due to a single ion channel. 



fractions of the bulk kinetic energy that go into in the electrons and the mag- 
netic field, ee and eg respectively, are not free and independent parameters 
of coUisionless shock theory. Most likely they represent interconnected parts 
of the same process. 

4. In the case of a weak or no upstream magnetic field, the Weibel- 
like two-stream instability is able to provide the necessary electromagnetic 
fields. We have shown here that the coUisionless shocked region is relatively 
thin, and we suggest that the non-thermal radiation observed from GRB 
afterglows and relativistic jets in general is emitted from such a relatively 
thin shell. 

It is clear that the non-thermal electron acceleration, the ion current fila- 
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mentation, the magnetic field amplification/generation, and hence the strong 
non-thermal radiation from the shock, is beyond the explanatory capacity of 
MHD. Whether or not the relativistic MHD jump conditions become valid 
on any larger scale is not possible to decide from the simulations presented 
in this chapter. 



4.4 Conclusions 76 



5. THE GLOBAL STRUCTURE OF COLLISIONLESS 

SHOCKS 



In collisionless shocks the mean free path of a particle is greater than the 
extent of the shock interface. Hence the particle distribution functions are 
highly anisotropic and one cannot study them using fluid methods. Rather, 
the dominant means of collision is indirect, mediated by the collective elec- 
tromagnetic field. In Chapter 3 it was demonstrated that in collisionless 
shocks, where no large scale magnetic field exists beforehand, the result- 
ing electromagnetic field is largely dictated by the evolution of two-stream 
instabilities. In this chapter I study global charged particle dynamics in 
relativistic collisionless e^e" pair-plasma shocks numerically, using three di- 
mensional computer experiments with up to ~ 10^ particles, and present 
results on the application and limitations of two dimensional simulations for 
the study of the global structure in ion-electron plasmas. The pair plasma 
simulations are advanced to a quasi-steady state, and compared to a fluid 
model. There is good agreement with the fluid model for physically rele- 
vant quantities, such as bulk density and momentum, which shows that the 
evolution can be extrapolated to larger timescalcs. We find that the two- 
stream instability decays over 50-100 electron skin depths and hence for a 
pair plasma shock remains firmly in the microphysical domain. This type of 
microphysical experiment may be used to determine empirically an effective 
equation of state in a collisionless shock, and the necessary sink and source 
terms that describe the conversion of kinetic to magnetic energy due to the 
two-stream instability, which could then be implemented in global fiuid mod- 
els, leading to more accurate large scale simulations of phenomena such as 
gamma ray bursts and relativistic jets from AGN's, where coUisionless shocks 
are of importance. The technique would be similar in spirit to the role played 
by subgrid models in understanding large scale turbulence, where models of 
the small scale behaviour are integrated into the fluid simulations to extend 
the dynamical range. 
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5.1 Introduction 

Three dimensional particle-in-cell experiments of ion-electron coUisionless 
shocks with open boundaries in the streaming direction have been considered 
in Chapters 3 & 4, and by Fredriksen et al. [25, 26], Hededal et al, [32, 33] and 
Nishikawa et al. [60], but as shown in Chapter 4 the estimated true extent 
of a coUisionless ion-electron shock is much larger than the computational 
domains that have been considered up to now. In Chapter 4 I found that 
for a mass ratio of mj/mg = 16 the shock extends at least 1500 ion skin 
depths. Unfortunately, with current computer technology there is no hope 
of performing three dimensional experiments that resolve scales all the way 
from sub electron-skin depths to 1500 ion skin depths. 

In Chapter 3 a qualitative difference between ion-electron shocks and pair 
plasma shocks was noted: In the case of an ion-electron plasma the heavier 
ions disrupt the electron channels and the electrons form a Dcbyc shielding 
cloud around the ion channels. This cloud of electrons stabilises the ion 
channels; indeed this is why the channels survive significantly longer and ions 
from the upstream and downstream mediums interact less. The consequence 
is that thermalisation of the ions and decay of the two-stream instability in 
an ion-electron dominated shock interface happens on a fundamentally longer 
time scale than in a shock interface dominated by a pair plasma. 

Even though full three dimensional ab initio experiments of ion-electron 
shocks are out of reach, that is not so for pair plasma shocks. In a pair plasma 
the electrons and positrons generate channels on the same time scale, and 
with no shielding they are quickly disrupted. In terms of the electron skin 
depth time scale Wp/c, thermalisation is faster than in the ion-electron case. 
Furthermore, the electrons and positrons have the same mass, and therefore 
many more skin depths can be resolved in a single box. 

5.2 CoUisionless Pair Plasma Shocks 

The two-stream instability deflects particles in the transverse direction to 
the flow, and to correctly describe a coUisionless shock the model has to 
be fully multi-dimensional. It was shown in Chapter 3, that the current 
channels merge in a self similar process generating a powerlaw distributed 
magnetic fleld. Medvcdcv ct al. [52] investigated the problem theoretically 
and demonstrated that the magnetic fleld correlation length in the shock 
interface grows with the speed of light for relativistic shocks. 

It is therefore necessary that our computational domain perpendicular to 
the shock is comparable in size to the longest two-stream unstable regions in 
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Fig. 5.1: The evolution of the total magnetic energy in the box as a function of 
time for the three runs. 



the box. Otherwise the process may be hmited numerically, and with periodic 
boundaries perpendicular to the shock interface the experiment reaches a 
state containing just a few self interacting current channels; then it cannot 
be entirely clear if the saturation of the magnetic field, decay of the current 
channels, and thermalisation of the particles happen due to numerical or 
physical effects. 

The basic setup of the numerical experiment has already been described 
in the previous chapters. In this section I consider three variants of the 
same experiment. All are pure pair plasmas. Initially the dense downstream 
medium is at rest. The density jump is a factor 3, and the inflow Lorentz 
factor of the upstream medium is 3. The only differences between the setups 
are in the box sizes and the plasma frequencies considered. They all contain 
initially 8 particles per species per cell in the medium at rest. In the main 
experiment. A, the box size is nx x ny x nz = 80 x 80 x 2800 and the plasma 
frequency is Up = 0.42. In the two complementary experiments, B and C, 
the box sizes are 160 x 160 x 1400 and 80 x 80 x 2800, with plasma frequencies 
of Up = 0.21 and Up = 0.105 respectively. 
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In experiment A plasma oscillations are resolved with only 2.4 cells, which 
is close to the Nyquist frequency of the grid, and indeed there is a higher 
level of small scale numerical noise than in experiments B and C. While 
it would have been preferable with an experiment with a smaller plasma 
oscillation frequency, the presented runs are at the limit of current computer 
capacities, containing up to a billion particles, and only experiment A settles 
to a steady state with selfsimilar evolution. Experiment B and C are used to 
vahdate the behaviour for early times. In fact the first stages of a thermalised 
shock is observed in experiment B, but it does not separate into different 
states before reaching the edge of the box. Specifically, the evolutions in 
averaged current and mass densities are equal for early times in the different 
experiments. Figure 5.1 shows the evolution of the total magnetic energy in 
the box. There is a clear difference in the initial level of fiuctuations between 
the experiments, due to the difference in plasma oscillation frequencies, but 
the growth rate is the same for the three cases, and the experiments show the 
same late time behaviour. We can separate the evolution in three phases: 
First an initial inflow phase, where the particles have not yet undergone 
the two-stream instability. At around t — 10a;~^ the two-stream instability 
commences, and it is saturated at around t — 40 — 100 o;"^. From then on 
the region containing shocked material and a diffuse turbulent magnetic field 
expands, while a reverse and a forward shock are formed, and a slow rise in 
the total magnetic field is seen. Both at the forward and the reverse shock 
interface a permanent two-stream instability is observed. 

The large scale structure at t — llOOo;"^ for experiment A is shown 
in Fig. 5.2. The average density profile has a pile-up of matter from Z = 
350 — 700, this is the shocked area. In Fig. 5.2 (a) profiles for earlier and 
later times have been overplotted to illustrate how the shocked area expands 
with time, with the forward shock to the right moving faster than the reverse 
shock to the left. Panel (b) in Fig. 5.2 shows the relative amount of power 
in the strongest mode of a Fourier transform of the transverse magnetic 
field compared to the integrated power. It has been calculated by taking 
the Fourier transform of the magnetic field in the x — y plane, finding the 
largest amplitude mode, and then dividing that with the integrated power. 
In the two-streaming shock interface, the largest transverse scale of the box 
is dominant, and the ratio is close to 1, while in the centre of the shocked 
medium the magnetic field has decayed and power is distributed over a range 
of scales. In (c) the magnetic energy in the field is plotted and the two 
shock interfaces are clearly seen to be separated. From Fig. 5. 2 (c) it can be 
estimated that the two-stream unstable regions have a width of between 50 
and 100 electron skin depths. 

In order to validate that the jump conditions are satisfied I have used the 
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Fig. 5.2: The large scale structure at t = 1100 oj^^. (a) The average density of elec- 
trons in an x—y slice, with similar plots for t = 1016 o;"^ and t = 1184 o;"^ 
as dashed lines, (b) The relative amount of power in the dominating mode 
compared to the total power in a two dimensional Fourier transform of 
the transverse magnetic field, (c) Average energy density in the magnetic 
field, (d) Average four velocity in the z-direction and a similar profile 
from an HD run. (e) Density from an HD run. 



fluid code described in Chapter 2 to setup a similar shock problem, but in 
fluid dynamics. I have chosen a relativistic equation of state with F = 4/3 
and a density jump of 3, an inflow four velocity of 1.7 and a uniform pressure 
of P = 0.2. The velocities and densities are taken in accordance with the 
unperturbed states seen in experiment A around t — 1100 o;"^. A priori 
it is not clear how to measure the pressure, since there are contributions 
from the random magnetic field, the two-stream generated magnetic field, 
heating from backscattered particles to the left of the shock, and heating 
from particles that were not scattered to the right of the shock. Instead of 
trying to measure an ill defined pressure, I have chosen the pressure P such 
that the maximum density in the right shock wave is in accordance with the 
particle data. This is a reasonable approach, because by fixing the pressure 
according to one single parameter, we see good agreement between the two 
models for the velocity profile and density profile in the other parts of the 
shock wave. Furthermore the shock profile is seen to move with a velocity of 
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0.6 (Approximately 50 skin depths in 84 skin depth times, see Fig. 5.2 (a)), 
while the corresponding fluid shock, in good agreement, is moving with a 
velocity of 0.65. Naturally, the profiles and bulk velocities do not correspond 
exactly, since not all facets of the particle experiment are reflected in the fluid 
shock. The hydrodynamical experiment does not include magnetic fields in 
any way, since no magnetic fields are present in the initial state. In contrast to 
that, in experiment A, strong magnetic field generation at the discontinuities 
in the velocity profile is seen. Moreover, the discontinuities are still rather 
smooth in experiment A, due to the collisionless nature of the shock. To 
take an example, some of the upstream particles, coming from the left, do 
not scatter in the shocked region. As a result there is a smooth transition 
in velocity and density at the forward shock front to the right. If the box 
size was larger, and the shock could be followed for longer times, making the 
extent of the shocked medium larger compared to the two-stream interfaces 
and effective mean free path, this smooth transition would stay constant and 
the solution would converge even better to that of a fiuid shock. A detailed 
analysis between larger experiments running for longer timescales and MHD 
shock tubes with a range of magnetic field configurations has to be done to 
fully understand the implications for the jump conditions of the two-stream 
instability; but this first experiment has shown that indeed it is possible to 
recreate the fluid representation ab initio of a pair plasma collisionless shock 
using a particle code and working from first principles. 

Recently similar experiments were presented by Spitkovsky [74], and his 
results are in agreement with our findings. 

The relatively short thermalisation length observed in this experiment 
implies that in the case of a pair plasma shock for astrophysical applica- 
tions, the two-stream instability remains a purely microphysical phenomenon, 
which probably has little impact on any observed quantities in astrophysical 
shocks, simply due to the small volume in which it takes place. 

5.3 Collisionless Ion-Electron Shocks and Limits to Current 

Experiments 

The computer experiments presented in last section are only a few of a series 
of large scale experiments that I have performed during the last year. The aim 
is to understand the global structure in ion-electron dominated collisionless 
shocks, by performing a series of three-dimensional experiments with lower 
and lower mass ratios in order to finally obtain a thermalised profile, and 
furthermore with this body of experiments to be able to rescale the results 
to realistic mass ratios. But even using mass ratios down to mi/rrie = 4, with 
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current computer limitations of around a billion particles, it is not possible 
to reach a state in the experiments, where both the ions and electrons fully 
thermalise and two shock interfaces emerge. 

A related but possible project is instead to understand coUisionless shocks 
in two dimensions. Two-dimensional collisionless shocks have been consid- 
ered for some time in the literature (e.g. Califano et al. [14] and Kazimura et 
al. [38]), but until now no large scale 2D simulations have been made with 
open boundaries (though see Medvedev et al. [52] for an experiment with 
periodic boundaries in the streaming direction, and the work by Hededal 
[31]). 

To compare 2D and 3D simulations, in this section I present two ex- 
periments with exactly the same initial conditions as the 3D experiment 
considered in Chapter 4 (from now on experiment C). They both have an 
inflow Lorentz boost of F = 15, an electron plasma frequency of cUp^g — 0.3 
and a mass ratio of rrii/me — 16. I use 8 particles per species per cell. 
Experiment A has the same box size transverse to the flow as experiment 
C with nx X nz = 125 x 4000, while experiment B is much wider with 
nxxnz^ 1600 x 4000. 

I have selected the two experiments to address the following questions: 1) 
How much do 2D and 3D experiments differ, both quantitatively but also in 
the underlying physical process. 2) What is the impact of the narrow boxes 
that have been considered until now. In nature a collisionless shock is much 
wider than the shock interface is long, and during the instability, the different 
regions, by causality, can only interact with a finite area of the shock front. 
But in some of the 3D experiments, to grasp the streaming nature of the 
shock properly, the boxes have been far too small in the transverse direction 
of the shock. 

In Fig. 5.3 the size of the current densities for the three experiments are 
shown at time t — 125 a; ~g. We see basic agreement between experiment A 
and experiment C for the morphology, thoTigh the ion current channels in 
experiment C arc thicker and more well structured. 

Comparing experiment A and B we see that the larger box size leads 
to a much more complex picture than the simple idea of current channels 
streaming strictly in the direction of the flow. There is a dazzling array of 
interactions going on, with nontrivial interactions between the channels. In 
the lower right of experiment B there is an almost square formed complex 
of current channels, which in itself is wider than both experiments A and 
C. The filamentary structure is sustained, but in contrast to the simple toy 
model presented in Chapter 3 and by Medvedev et al. [52] , one cannot speak 
of a merging hierarchy of ion channels, and the lifetime of an individual 
channel is quite small. 
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Fig. 5.3: From top to bottom: The current density of the ions and electrons of 
experiment B, experiment A and the averaged current density in the y- 
direction of experiment C, reported on in Chapter 4. The dashed Unes 
indicate the region used for constructing particle distribution functions. 
The figure is shown with the correct aspect ratio and the snapshots are 
taken at t = 125 w^g. Length units are given in electron skin depths. 



The process is initially ignited by the two-stream instability, and in exper- 
iment B abundant examples of forming channels may be found, not only at 
the initial interface at Z = 200, but also downstream of the first generation, 
due to two-stream like configurations of the magnetic field. Moreover we 
observe direct merging and head on collisions between the individual chan- 
nels. Counteracting this process and partly responsible for the decay of ion 
channels is the electric potential. Near the centre of the channel there is a 
strong over density of positive charges and even after taking into account a 
relativistic time dilation factor - or equivalently - the self generated magnetic 
fields of the channels, they will "explode". In two dimensions they leave a 
cone-like structure containing two trails of ions, with some symmetry, since 
everything happens at the speed of light, while in three dimensions a ring- 
like structure is created. A three-dimensional version of this explosion may 
be seen in Fig. 4. IB, where, except for the helix structure, everything is sta- 
bilised and symmetric due to selfinteraction of the channel. This electrostatic 
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Fig. 5.4: Particle distribution function for the electrons in a slice around Z = 
400 u;~g . To the left is shown the PDF for experiment B and to the right 
the PDF for experiment C. The PDF for experiment A is identical to B. 



explosion makes the evolution of collisionless shocks even more dynamic and 
intermittent. It is important to point out that the effect depends on the effec- 
tiveness by which the electrons Debye shield the ion channels, and therefore 
on the mass ratio of the experiment. At higher, more realistic mass ratios, 
the shielding is more effective and the timescale for the breakup of the ion 
channel longer. 

The relatively short time scale of experiment B does not make it possible 
to assess what the long time implications are of this richer structure, but 
does show that it is important to have an adequate resolution transverse to 
the flow, and not only in the streaming direction, if a full understanding 
of collisionless ion-electron shocks is to be obtained. The highly dynamic 
nature and evolution along the streaming direction also show that streaming 
experiments with open boundaries are essential to understand the state of 
the plasma far downstream of a collisionless shock interface. 

It is important to understand if there are differences in the morphology of 
the currents observed in experiments with two- and three-dimensional shocks. 
But in order to make a more formal and quantitative investigation we have to 
look at the particle statistics. I have measured the particle distribution func- 
tions (PDFs) for the electrons in a slice of the domain delimited in Fig. 5.3 
by dashed hnes for the three experiments. The shce has been selected to lie 
at the point in the shock where the electrons are on the brink of merging 
to one single continuous population, but still the form of the PDF is domi- 
nated by remnants of the initial upstream and downstream populations. The 
slope of the PDF indicated in the figure depends on the amount of heating 
in the populations. A warmer upstream population will be broader in phase 
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space, and consequently the maximum is lower, giving rise to a steeper slope. 
It should be emphasised that the perceived power law seen in the figure is 
but a consequence of the merging populations. It can easily be verified by 
noting that the "powerlaw" breaks at around — 15, the velocity of the 
instreaming population. We find perfect agreement between the PDFs of the 
two-dimensional experiment A and B, and they both have a slope index of 
2.1, while experiment C has a slope index of 1.55 (see Fig. 5.4). The reason 
for this significant difference in the heating rates of the electrons in two- 
dimensional shocks compared to three dimensional is the same as the reason 
for the lower lifetime of the ion channels in the two-dimensional experiments. 

Essentially we can understand the differences between two- and three- 
dimensional simulations by considering the real space available in the two 
cases. 1) There is a difference in the dynamics of the ion channels: If we make 
a transverse cut through the flow in the three-dimensional case, the channels 
may be likened to ID particles in a plane (see Fig. 3.1 in Chapter 3 for an 
illustration of this). We have carefully studied the time evolution of this ID 
gas of current channels in the experiment considered in Chapter 3. It is found 
that, if channels collide head on, they will generally not coalesce, and instead 
in many cases destabilise, because of the inertia of each channel, while in cases 
where the impact parameter is very low, or the two channels slightly miss each 
other, a much smoother collision process is initiated by the in-spiralling of 
the two channels, ultimately leading to the formation of a single channel. In 
a two-dimensional simulation, though, the two-dimensional transverse plane 
reduces to a one dimensional line, and consequently two merging channels 
will always collide head on, making coalescence more difficult, the transverse 
velocities, i.e. the temperature, of the ions higher, and the lifetime of the 
channels lower. 2) There is a difference in the dynamics of the electrons: The 
electrons Debye shield the ion channels, and move generally in accordance 
with the toy model described in Chapter 4. In the three-dimensional case the 
paths of two such electrons have been depicted in Fig. 4.1 in Chapter 4. In 
general the electrons do not move exactly through the centre of the current 
channel, but instead traverse some kind of complicated ellipsoidal path. In 
the two-dimensional case though, because there is no "third dimension" to 
miss the centre in, the electrons have to go right through the centre of the ion 
channel, and gain the maximum amount of acceleration. The acceleration 
is to a large extend potential, and the electron looses most of the energy 
climbing out of the potential well, but because of the time dependence of 
the fields, statistically there will be some momentum transfer. In the two- 
dimensional case the electrons have to pass through the local minimum of 
the electrostatic potential, and hence the heating is more effective than in 
the three-dimensional case. 
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The experiments presented in this section have shown that, while the ba- 
sic morphology and dynamics do carry over from three to two dimensions, 
there are some quantitative differences. The heating of electrons and ions is 
more effective and the two-stream generated ion channels are not as stable in 
two as in three-dimensional experiments. Nonetheless, if we take these dif- 
ferences into consideration when interpreting two-dimensional experiments, 
these experiments are still the most promising tools for understanding the 
global structure of ion-electron coUisionless shocks. Prom the above discus- 
sion it is clear that the extent of the two-streaming region in an ion-electron 
coUisionless shock, as inferred from future two-dimensional experiments of 
the global shock structure, will most likely be smaller than in the case of a 
three-dimensional shock. This conclusion may be drawn directly from the 
higher heating rate alone, as observed in Fig. 5.4. 

5.4 Conclusions 

CoUisionless shocks arise in many astrophysical objects and the correct \in- 
derstanding of relativistic coUisionless shocks have implications for the obser- 
vations of outUows from compact objects, such as gamma ray burst afterglows 
and relativistic jets. We have seen in the preceding chapters that magnetic 
field generation and particle acceleration are integral parts of coUisionless 
shocks in the case of weak or absent large scale magnetic fields. To under- 
stand the impact on observations it is essential to investigate how far down 
stream of the initial shock the two-stream unstable region extends. With 
this in mind I have, in the current chapter, discussed the global structure of 
coUisionless shocks. 

In the first part of the chapter I have presented a fully three-dimensional 
model of colliding pair plasmas using a particle-in-cell code, and observed the 
thermalisation of the plasma due to the collective electromagnetic field, and 
the formation of a macro physical shock structure. Comparing the results to 
a fiuid simulation, with the same initial conditions, good agreement is found, 
implying that the full structure of the shock has been resolved. 

Crucially, I have estimated that the decay of the two-streaming region and 
subsequent thermalisation happens over 50-100 electron skin depths. Thus, 
for the specific case we have considered it renders the two-stream instability 
for pair plasmas completely microphysical. Hence, the two-stream instability 
in coUisionless shocks comprised purely of leptonic matter may have few 
direct observational consequences. 

In the second part of the chapter I have considered the global structure of 
ion-electron dominated coUisionless shocks. With current computer capaci- 
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ties it is impossible to correctly model the global structure of an ion-electron 
shock in three dimensions. Two-dimensional coUisionless shocks remain a 
promising alternative, and I have investigated their applicability in under- 
standing three-dimensional models. It has been shown that while indeed 
two-dimensional shocks, for the time being, are our best hope to grasp nu- 
merically the global structure of ion-electron coUisionless shocks, there are 
some differences, and caution should be voiced in directly generalising results 
from two-dimensional experiments to three dimensions. The ion channels 
that form due to the two-stream instability are less stable, and the heating 
rate of the electrons is higher. Both factors contribute to a faster thermal- 
isation than what can be expected from three-dimensional experiments in 
the future, and hence give rise to an underestimation of the extent of the 
two-stream unstable region. Nonetheless, the overall physical picture is the 
same, and these differences can be taken into account. 



6. A NEXT GENERATION PIC CODE 



6.1 Introduction 

Over the last couple of years the Copenhagen group has been using PIC mod- 
els that include electromagnetic fields and charged particles to understand the 
plasma microphysics of coUisionless shocks [25, 26, 32, 33]. It has turned out 
to be a very successful tool, but it is still limited in the scope of phenomena 
that may be addressed. Even though a large class of astrophysical environ- 
ments are indeed coUisionless, scattering and collision processes do play an 
important role in several key scenarios. Examples are given below. Another 
key ingredient, which has been missing in charged particle simulations, is 
a full treatment of photon propagation. It can be argued that photons are 
represented directly on the mesh by electromagnetic waves, which certainly 
is correct. But the mesh can only represent waves with frequencies smaller 
than the Nyquist frequency. The physical length of a typical cell has in our 
apphcations typically been 10^ — 10^ cm and hence it is clear that only low 
frequency radio waves can be represented. High frequency photons have to 
be implemented as particles that propagate through the box and interact, 
either indirectly through messenger fields on the mesh, or directly with other 
particles. A valuable consequence of modeling the detailed photon trans- 
port is that extraction of electromagnetic spectra is trivial. Even in cases 
where the photon field is only a passive participant, this fact should not be 
underestimated as it enables direct comparison with observations. 

There exists Monte Carlo based particle codes (see e.g. [76] and references 
therein) that address various particle interactions, but one of their main 
shortcomings is the poor spatial resolution. This makes it impossible to 
couple the particle aspects to a self consistent evolution of the plasma. 

Our goal has been to develop a framework where both electromagnetic 
fields and scattering processes arc included in a consistent way. We can 
then correctly model the plasma physics and the radiative dynamics. The 
scattering processes include, but are not limited to, simple particle-particle 
scattering, decay and annihilation/creation processes. Our new code is not 
limited in any way to charged particles, but can also include neutrals such 
as photons and neutrons. 
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In the next section we describe some of the physics that can be addressed 
with this new code. In section 6.2 we discuss how the code has been imple- 
mented, the general framework and in detail which physical processes that 
are currently implemented. In section 6.3 we present the first results in the 
form of a toy experiment that we have performed to validate the code. In 
the last section 6.4 we summarize. 

6.1.1 Motivation 

Before we continue and describe in detail the methods, physics and test 
problems we have implemented and used, it is important to consider the 
general class of scenarios we have had in mind as motivation for developing 
the code. There are several key objects, where only the bulk dynamics is 
understood, and we are lacking detailed understanding of the microphysics. 

Internal shocks in gamma ray bursts 

In the internal/external GRB shock model, the burst of gamma-rays is be- 
lieved to be generated when relativistic shells collide and dissipate their rel- 
ative bulk energy [56, 70]. The nature of the radiation is presumably inverse 
Compton scattering and synchrotron radiation. Particle/photon interactions 
might also play an important role in the very early afterglow as suggested 
by [9, 77]: Even though the medium that surrounds the burst (ISM or wind) 
is optically very thin to gamma-rays, a tiny fraction of the gamma-rays will 
Compton scatter on the surrounding plasma particles. This opens up for 
the possibility of pair-creation between back scattered and outgoing gamma- 
rays. The creation of pairs may increase the rate of back scattered photons 
in a run- away process [75]. The Compton scattering may accelerate the 
pair-plasma through the surrounding medium with many complicated and 
non-linear effects, including streaming plasma instabilities and electromag- 
netic field generation. Hence, it is crucial that plasma simulations of internal 
GRB plasma shocks include lepton-photon interactions. 

Solar corona and the solar wind 

Space weather (defined as the interaction of the solar wind on the Earth) 
is in high focus for several reasons. Not only is the Sun our closest star, 
providing us with invaluable data for stellar modeling, but coronal mass 
ejections from the Sun potentially have impact on our every day life. The 
strong plasma outflows from the sun can induce large electrical discharges in 
the Earths ionosphere. This may disrupt the complex power grids on Earth, 
causing rolling blakcouts such as the one in Canada and North America 
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in 1989. Also high-energy particles can be hazardous to astronauts and 
airline passengers. Computer simulations have provided a successful way 
of obtaining insight in these complex plasma physical processes. However, 
in the solar coronal and in the solar wind plasma out to distances beyond 
the earth orbit, difficulties arise in finding the right formalism to describe 
the plasma. Neither a collisionless model based on the Vlasov equation nor 
an MHD fluid model provides a adequate framework for investigation. The 
problem has already been studied using three dimensional PIC simulations 
but without taking collisions into account (e.g. [13, 35]). 

The corona of compact objects 

The bulk dynamics of accreting compact objects have been modeled for many 
years using fluid based simulations (e.g. [5] and references therein). Never- 
theless, it has been a persistent problem to extract information about the 
radiating processes. Furthermore in the corona the MHD approximation 
becomes dubious, just as in the solar corona. The environment around a 
compact object is much more energetic than the solar corona, and therefore 
radiative scattering processes play an important role. Pair production is also 
believed to be abundant. Using our new code it would be possible to model 
a small sub box of the corona. The main problem here - as in most numer- 
ical implementations - is to come up with realistic boundaries for the local 
model. A shearing box approach may be appropriate, but in fact we can do 
even better. 

The size of a stellar mass black hole is around 10^ cm. In a fluid simulation 
we want to model the accretion disk-compact object system out to hundreds 
of radii of the compact object. The normal approach is to use a non uniform 
mesh. Nonetheless, the Courant criterion, which determines the time step, is 
still limited by the sound crossing time of the compact object. I.e. the time 
step is limited by the size of the innermost (and smallest) cells in the mesh. 
The very small time step corresponds to those found in a typical particle 
simulation, where the strict time step arises from the need to resolve plasma 
oscillations. Hence data from an MHD simulation could provide temporally 
well resolved fluxes on the boundaries of the much smaller sub box containing 
the particle simulation. 

In this sense the particle simulation will act as a probe or thermometer 
of the fluid model. The particle model includes the full microphysics in a 
realistic manner and most importantly includes photon transport. Realistic 
spectra could be obtained indirectly from the fluid model, testing fluid theory 
against observations. We have already worked on interfacing fluid models 
with the old PIC code. 
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Pre-acceleration in cosmic ray acceleration 

Accepting Fermi acceleration as a viable mechanism for accelerating electrons 
and creating the the non-thermal cosmic ray spectrum is still left with some 
big unanswered questions. One is that the Fermi mechanism requires injec- 
tion of high energy electrons while still keeping a large, low-energy population 
to sustain the magnetic turbulence. Hence, a pre-acceleration mechanism 
needs to be explained. 

The shocks in supernova remnants are believed to be cosmic ray accelera- 
tors. However, the Fermi acceleration process in shocks is still not understood 
from first principles but rely on assumptions on the electromagnetic scatter- 
ing mechanism. PIC codes would seem ideal in exploring the mechanism 
from first principles, since they include field generation mechanisms and the 
back-reaction that the high energy particles have on this scattering agent. In 
Supernova remnants however, the mean free path for Coulomb collisions are 
comparable to the system and particle-particle interactions cannot be fully 
neglected. 

6.2 Implementation 

Implementing any state-of-the-art large scale numerical code is a big under- 
taking, and can easily end up taking several man years. We estimate that the 
final version of the next generation code will contain more than 50.000 lines 
of code. Starting in February this year, it has taken us three man months 
to implement the current incarnation of the code which has already grown 
to approximately 10.000 lines. Besides T. Haugb0lle and C. B. Hededal, the 
development is done together with A. Nordlund and J. T. Frederiksen. For- 
tunately we have a good tradition and expertise for numerical astrophysics 
in Copenhagen and we have been able to port different technical concepts 
and solutions from our suite of fluid codes and to a lesser extent from the 
old PIC code. The aim is to build an extremely scalable code that is able to 
run on thousands of CPUs on modern cluster architectures and utilize MPI 
as the inter node communication protocol. In this chapter we will not go 
further into technical details. Instead we will put emphasis on the important 
concepts and physics and how we have implemented these. 

6.2.1 Concepts 

The two fundamental objects in a particle- in-cell code are the mesh and the 
particles. Wc have adopted the solver and interpolation routines from the 
old PIC code to solve the Maxwell equations and find fiuxes and densities 
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on the mesh. The mesh is used to distribute messenger fields - such as the 
electromagnetic fields - and to calculate volume averaged fiuxes and densities 
of the particles The latter are used as source terms in the evolution of the 
messenger fields. The particles really represent an ensemble of particles and 
are often referred to as pseudoparticles [10] or large particles. A so called 
smoothing kernel describes the density distribution of a single pseudoparticle 
on the mesh. In our implementation the volume of a particle is comparable 
to a cell in the mesh. 

Pseudoparticles with variable weights 

The concept of pseudoparticles is introduced since the "real space" particle 
density easily exceeds any number that is computationally reasonable (i.e. 
of the order of a billion particles). The pseudoparticle charge to mass ratio 
is kept the same as the ratio for a single particle. 

In ordinary PIC codes the weight of each pseudoparticle of a given species 
is kept constant throughout the simulation. The benefit is a simple code and 
a unique identity for each particle. The first is a convenience in the practi- 
cal implementation, the second important when understanding the detailed 
dynamics and history of a single particle. 

Notwithstanding possible conveniences, as detailed below in section 6.2.1, 
we have decided to improve this concept to a more dynamical implementa- 
tion where each pseudoparticle carries a individual weight. Particles are then 
allowed to merge and split up when a cell contains too many/few particles, 
or when particles are scattered. The concept is sometimes used in smooth 
particle hydrodynamics (SPH), where different techniques have been pro- 
posed for the splitting and merging of particles. It is both used to adjust 
the density of individual particles [11] and in the conversion of gas- to star 
particles in galaxy formation models [28]. An important quality of SPH is its 
adaptive resolution capabilities. These are important in the description of 
collapsing self gravitating systems, ranging from core collapse supernovae to 
the formation of galaxy clusters, scenarios where matter is collapsing many 
orders of magnitude, and therefore the smoothing length or volume of the 
individual particles is readjusted accordingly. Consequently, when splitting 
particles or adjusting the weights in an SPH code, it is important to match 
precisely the spatial density distribution of the parent particle to the spatial 
distribution of the child particles. In PIC codes, though, the spatial size or 
smoothing parameter of an individual particle is determined beforehand by 
the mesh spacing. This is reasonable since we are not interested in adaptive 
resolution but rather a kinetic description of the plasma dynamics. Split- 
ting a parent particle with weight Wp into child particles with weights is 
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therefor trivial. The requirements of conservation of mass and four velocity 
together with conservation of the density and flux distribution in the box, 
can all be satisfied by setting 

n 

Wp = ^wl, ep^el, -fpVp = -flvl , (6.1) 

i=l 

since the smoothing kernel is determined by the mesh spacing, not the mass 
of the individual particle. 

The merging or renormalization of pseudoparticles requires a much more 
thorough analysis. Up to now we have investigated two schemes, one that 
respects conservation of mass, energy and four velocity by merging three 
particles into two at a time, and one where only mass, energy and average 
direction is conserved by merging two particles into one particle. While 
these schemes probably work well for approximately thermal distributions, 
it will easily give rise to a large numerical heating when considering head on 
beam collisions. We believe it can be improved by first selecting a random 
"merger particle" and then find other particles in the local cell, that are 
close to the merger particle in momentum space. A more radical approach 
is to resample the full phase distribution in a cell every time the number 
density becomes above a certain threshold. Nevertheless, it requires testing 
of different extreme situations to find the optimal method to merge particles, 
and it is still a work in progress. 

To obtain the results that we present in section 6.3, we ran the code 
without merging of the pseudoparticles activated. 

Scattering processes and splitting of particles 

In Monte Carlo based particle codes the generic way to compute an interac- 
tion is first to calculate the probability for the interaction Ps, then compute 
a random number a. If a < Ps then the full pseudoparticle is scattered oth- 
erwise nothing happens. This probabilistic approach is numerically rather 
efficient and simple to implement, but it can be noisy, especially when very 
few particles are present in a cell. In large particle Monte Carlo codes the 
typical cell contains up to 10"^ particles per species per cell (hence "large par- 
ticle"). In our PIC code typical numbers are 10^ — 10^ particles per species 
per cell, since we need many cells to resolve the plasma dynamics. For our 
requirements the probabilistic approach would result in an unacceptable level 
of noise. For example, in a beam experiment the spectra of the first gener- 
ation of scattered particles may come out relatively precise, but the spectra 
of higher generation scattered particles (i.e. particles that are scattered more 
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Fig. 6.1: To implement the scattering of two pseudoparticles we transform to the 
rest frame of the target particle (shown as red/light gray) and computes 
the probability P{n) that a single incident particle (shown as blue/dark 
gray) during a time step At is scattered on the n target particles. If the 
incident particle has weight m, then k = P{n)m, particles will interact 
and two new pseudoparticles are created. 



than once) will come out with poor resolution or require an excessive amount 
of particles. Another well known consequence of the probabilistic approach 
is that for a given experiment the precision goes in the best case inversely 
proportional to the square root of the number of particles used in the exper- 
iment. To increase effective spectral resolution we have instead decided to 
take a more direct approach. For simplicity we will here describe the method 
for a two-particle interaction, and disregard all factors converting code units 
to physical units. For example, the weight of a pseudoparticle is proportional 
to the number of physical particles in the pseudoparticle. Although, these 
prefactors all represent trivial conversions of units, they must be taken into 
account in the real code. 

Consider a single cell containing a single pseudoparticle (red) with weight 
Wt = n and a single pseudoparticle (purple) with weight Wi = m, where 
n > m (see Fig. 6.1). We first select the red particle as the target, since 
n > m, and the purple as the incident particle. We then transform the 
four velocity of the incident particle to the rest frame of the target particle, 
and calculate the total cross section at of the interaction. Conceptually we 
consider the process as a single incident particle approaching a slab of the 
target particle. The number density of target particles in the slab can be 
calculated from the weight Wt as pt = Wt/AV, where AV = AxAyAz is 
the volume of a single cell. Given the number density the probability that a 
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single incident particle is scattered per unit length is 

Pi = PtCTt = ^ • (6.2) 

During a time step At the incident particle travels Al = ViAt, and the 
probability that a single incident particle is scattered then becomes 



P5 = l-exp[-P;AZ] 

WtatViAt 



1 — exp 



AV 



(6.3) 



The weight of the incident pseudoparticle is Wi = m. Pseudoparticles repre- 
sent an ensemble of particles. Therefore Ps is the fraction of incident particles 
that are scattered on the target. To model the process we create two new 
particles with weight Wnew = WiPs = k. Given the detailed interaction, we 
can calculate the theoretical angular distribution of scattered particles in 
accordance with the differential scattering cross section. Drawing from this 
distribution we find the momentum and energy of the new scattered particles. 
The weights of the target and incident particles are decreased to Wt — n — k 
and Wi = m — k respectively (see Fig. 6.1). 

Our method faithfully represents the actual physics even for small cross 
sections. However, if all the particles are allowed to interact, the number of 
particles in the box will increase at least proportionally to the total number of 
particles squared. This is potentially a computational run away. Normally 
we will have on the order of up to 100 particles per species per cell, but 
to be computationally efficient we only calculate interactions for a subset 
of the particles in a cell. This subset is chosen at random according to 
an arbitrary distribution we are free to select. If the probability that two 
particles are selected for scattering in a given time step is Q then the traveling 
length Al simply has to be adjusted as Al/Q. If this arbitrary distribution 
is chosen cleverly, the particles with the largest cross section arc actually 
the ones selected most often for scattering, and everything ends up as a 
balanced manner: We only calculate the full cross section and scattering 
as often as needed, and the computational load that is given to a certain 
particle is proportional to the probability of that particle to scatter. We 
rely on the merging of particles as described above to avoid the copious 
production of pseudoparticles. Every time the number of pseudoparticles in 
a given cell crosses a threshold, pseudoparticles are merged and this way the 
computational load per cell is kept within a given range. 
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6.2.2 Neutron decay 

Free neutrons not bound in a nucleus will decay with a half-life a little longer 
than ten minutes. The neutron decays into an electron and a proton and an 
electron antineutrino to satisfy lepton number conservation 

n — > p + e~ + Pg (6.4) 

The rest mass difference of the process (0.78 MeV) goes into kinetic energy 
of the proton, electron and neutrino. Let the neutron lifetime be r in code 
units. If r is comparable to or less than a typical time step, then practically 
all neutrons decay in one iteration, and it is irrelevant to include them. If 
T is much larger than the total runtime, the neutron can be considered a 
stable particle (unless the neutron density in the box is much larger than 
the proton- or electron density). If instead r ~ a At where a ~ 100, then 
we can select a fraction / of the pseudoparticle neutrons in each cell and let 
them decay. This is done in an analogous manner to the generic scattering 
process described above in section 6.2.1. The weight of the selected neutron 
is decreased with a factor 

(6.5) 

where 7 is the Lorentz boost of the neutron pseudoparticle and / is chosen to 
give reasonable values for the decrease in the weight. At the same time a pair 
of electron and proton pseudoparticles is created with the same weight. The 
generated particles share the excess mass of the process (where the neutrino 
is neglected for now, but could be included in the future). The momenta are 
selected to give an isotropic distribution in the rest frame of the decaying 
neutron. 



exp 



/At 

7r 



6.2.3 Compton scattering 

Here we briefly describe a specific physical scattering mechanism which has 
already been implemented in the code, namely Compton scattering. 

Compton scattering is the rclativistic generalization of the classical Thomp- 
son scattering process, where a low energy photon scatters of on a free elec- 
tron. In the rest frame of the electron, the photon changes direction and 
loses energy to the electron which is set in motion. The cross section for 
Thompson scattering is [72] 

(^T = > (6-6) 

where vq = e^/(mc^) is called the classical electron radius. The Thompson 
scattering approximation is valid as long as the photon energy is much lower 



6.2 Implementation 



98 




Fig. 6.2: Schematic view of the Compton scattering process. Impinging on the 
electron, an incoming photon with energy ei is scattered into the angle 9 
with energy €2- In the initial rest-frame of the electron, the electron will 
be recoiled to conserve energy and momentum. 



than the electron rest mass hi/ <S mgC^ and the scattering can be regarded 
as elastic. For photon energies comparable to, or larger than, the electron 
rest mass, recoil effects must be taken into account. Measured in the electron 
rest frame we define ei as the photon energy before the scattering, 62 as the 
photon energy after the scattering and 9 the photon scattering angle (6.2). 
By conservation of energy and momentum one can show (e.g. [72]) that 

^ l + ^(l-cos^) ■ ^^-^^ 

The differential cross section as a function of scattering angle is given by the 
Klein-Nishina formula [34, 39] 

doc |4f!i + !^_sin^^V (6.8) 



The Klein-Nishina formula takes into account the relative intensity of 
scattered radiation, it incorporates the recoil factor, (or radiation pressure) 
and corrects for relativistic quantum mechanics. The total cross section is 
then 



3 



- ln(l + 2a;) y + — ln(l + 2x) - 



\ \^2x ' ' \ 2x ^ ' (l + 2a;)2 

(6.9) 

where x = hv/ {mc^). 
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Fig. 6.3: 3D scatter plot of a photon beam (black) passing through a cold pair 

plasma (gray). Left panel show initial setup where a photon beam is 
injected in the upward direction. Right panel shows how photons are 
scattered on the electron-positron pairs 



6.3 Results 

To test the new code and it capabilities in regard to the inclusion of colli- 
sions, we have implemented and tested a simple scenario involving Compton 
scattering. 

In the test setup, we place a thin layer of cold electron-positron pair 
plasma in the computational box. From the boundary, we inject a monochro- 
matic beam of photons all traveling perpendicular to the pair-layer (Fig. 6.3 
left panel). As the beam passes through the plasma layer, photons are scat- 
tered (Fig. 6.3 left panel). 

For each scattered photon we sample the weight of the photon and its 
direction (remembering that all particles arc pscudoparticles that represent 
whole groups of particles). Fig. 6.4 shows the theoretical cross section as 
function of scattering angle compared with the result from the simulations. 
Four plots for different energies of the incoming photon beam are shown. We 
find excellent agreement between the simulation results and the theoretical 
predictions. 
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Fig. 6.4: The theoretical Compton scattering differential cross section. We have 
performed a test experiment with an incoming laser beam on a very 
cold electron population. Over plotted the differential distribution is the 
theoretical curve according to Eqs. (6.8) and (6.7). 



6.4 Discussion 

A next generation PIC code that includes many different kinds of scattering 
process is under development. It will enable us to target problems that 
resides in the grey zone between the MHD and coUisionless plasma domains. 
This domain covers many astrophysical scenarios of great interest counting 
internal shocks in gamma-ray bursts, solar flares and magnetic substorms, 
compact relativistic objects, supernova remnants and many more. 

The concept of splitting/ merging particles and keeping individual weights 
of each particle carry many important features. Variable weights represent 
the true statistics of a scattering process in an optimal way compared to the 
Monte Carlo approach. Also, for MPI-parallelization it is crucial that the 
number of particles per cell is kept more or less constant to ensure an optimal 
CPU load-balancing. To localize calculations we are employing a sorting 
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algorithm that maintains neighboring particles on the mesh as neighbors 
in memory. This is not only good for parallelization, but also makes all 
computations very cache efficient; a crucial requirement on modern computer 
architectures. 

To test the infrastructure of the new code we have implemented Comp- 
ton scattering as a simple scattering mechanism. The first results are very 
promising in form of excellent agreement with the theoretical prediction. We 
note that a recent paper by [57] provide an interesting test suite for various 
kind of particle-photon interactions that can be tested in the future. Merging 
particles has not been satisfactorily implemented yet. Parallelization of code 
is still not there yet, and is necessary to obtain the capability of performing 
truly large scale experiments. In summary: Work has still to be done before 
we can start to investigate non trivial astrophysical scenarios, nevertheless 
solid progress has already been made 

This chapter has been written jointly by Christian Hededal and Troels 
Haugb0lle, reflecting the fact that the development process of the next gen- 
eration PIC code has been highly team based. Essentially everybody have 
contributed time and effort to every single source file of the code. It would 
not make sense to write the chapter separately, essentially repeating each 
other and reusing the same figures. 
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7. SUMMARY & CONCLUSIONS 



In the past chapters of this thesis I have presented different numerical meth- 
ods, as well as applications of the methods to a number of current problems 
in relativistic astrophysics. The thesis is logically structured into three parts, 
and below I would like to summarise the most important points of the work 
presented. 

In the first part (Chapter 2) I have presented the theoretical foundation 
and numerical implementation of a new general relativistic magnetohydro- 
dynamics code. 1 have derived a new form of the equations of motion, with 
global coordinates evolving the dynamical variables from the point of view 
of a local observer. When deriving the equations of motion, I have not made 
any assumptions about the background metric, so the design is ready to be 
coupled with methods solving the full Einstein equations. The code has been 
tested on a variety of demanding problems, and it has been demonstrated 
that it is able to deal with huge pressure and density gradients. The com- 
puter code is fully three-dimensional and parallelised and shows a spectacular 
performance on modern computer architectures exploiting up to 30% of the 
theoretical peak performance. It has been tested and verified to scale to 
hundreds of CPUs, making it possible to exploit massive supercomputers at 
national centres to the full extent. 

In the second part of the thesis (Chapters 3-5) I have presented important 
results in the understanding of coUisionless shocks using a charged relativis- 
tic particle-in-cell code. Together with Jacob Trier Prederiksen, Christian 
Hededal and Ake Nordlund I have investigated the fundamental consequences 
of the two-stream instability for observations of coUisionless shocks in gen- 
eral, and the implications for gamma ray afterglows in particular. In Chapter 
5 I extended our analysis and presented results on the global structure and 
transition of coUisionless shocks to fiuid shocks. 

In Chapter 3 we have shown that even in the absence of a magnetic 
field in the up-stream plasma, a small scale, fiuctuating, and predominantly 
transversal magnetic field is unavoidably generated by a two-stream insta- 
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bility reminiscent of the Weibel-instability. In the current experiments the 
magnetic energy density reaches a few percent of the energy density of the 
in-coming beam. 

In Chapter 4 we proposed an acceleration mechanism for electrons in 
ion-electron collisionless shocks. The acceleration mechanism is capable of 
creating a powerlaw electron distribution in a collisionless shocked region. 
The theoretical considerations were suggested by particle-in-cell computer 
experiments, which also allowed quantitative comparisons with the theoret- 
ical predictions. We have shown that the non-thermal acceleration of elec- 
trons is directly related to the ion current channels in the shock transition 
zone and is local in nature. The electrons are accelerated to a powerlaw 
in situ. Therefore the observed radiation field may be tied directly to the 
local conditions of the plasma and could be a strong handle on the physical 
processes. 

To understand the impact on observations it is essential to investigate 
how far down stream of the initial shock the two-stream unstable region 
extends. With this in mind I have analysed, in Chapter 5 the global struc- 
ture of collisionless shocks. I have presented three-dimensional experiments 
of colliding pair plasmas using the particle-in-cell code, and observed the 
thermalisation of the plasma, due to the collective electromagnetic field, and 
the formation of a macrophysical shock structure. Comparing the results 
to a fluid simulation, made using the code presented in Chapter 2, with 
the same initial conditions, good agreement is found, implying that the full 
structure of the shock has been resolved. I have estimated that the decay 
of the two-streaming region and subsequent thermalisation happens over 50- 
100 electron skin depths. Hence, the two-stream instability in collisionless 
shocks comprised purely of leptonic matter may have few direct observational 
consequences. 

In the second part of Chapter 5 I have considered the global structure 
of ion-electron dominated collisionless shocks. I have investigated the ap- 
plicability of global models using two-dimensional shocks - just possible 
with current computer technology - in the understanding of the complete 
three-dimensional shock structure It is demonstrated that caution should be 
observed in generalising results from two-dimensional experiments to three 
dimensions. In two dimensions the ion channels that form due to the two- 
stream instability are less stable, and the heating rate of the electrons is 
higher. Both factors contribute to a faster thermalisation than what may 
be expected from three-dimensional experiments in the future, and hence 
cause an underestimation of the extent of the two-stream unstable region. 
Nonetheless, the overall physical picture is the same, and these differences 
may be taken into account. 
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In the third part of the thesis (Chapter 6) together with Christian Hededal 
I have presented a new code under development by our group, which will en- 
able us to study not only charged particle dynamics, but also the propagation 
of neutral particles, such as photons and neutrons, as well as interactions be- 
tween these. 

The code is an extension of the current particle-in-cell code, and also 
solves the full Maxwell equations, but furthermore considers particle-particle 
interactions and microphysical processes, such as scattering, pair production, 
decay, and annihilation of particles. 

Especially the inclusion of photons and related radiative processes is im- 
portant. In the future we will be able to extract self consistent spectra from 
our numerical experiments, thereby gaining the ability to directly compare 
our models with observations. 

Even though the different tools presented in this thesis per se are not 
connected, they all revolve around the same physical problems. In Chapter 
5 we saw the first example of connecting the codes, to obtain different points 
of view on the same physical situation. In conclusion, and with a look to the 
future, I believe that the couphng of the GrMHD code with the new photon 
plasma code yields a great potential for obtaining realistic synthetic light 
curves from fluid simulations, connecting them directly with observations. 
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APPENDIX 



A. THE RELATIVISTIC MAXWELL DISTRIBUTION 



In this appendix I briefly consider the relativistic Maxwell distribution. When 
working with data from the particle code, we have often needed to assess 
how thermal or non thermal a given particle distribution function (PDF) for 
a subset of the particles is, and evaluate the temperature and the overall 
Lorentz boost of the population. Even if the particles are in fact thermally 
distributed, they can still be moving with an overall velocity u. To find 
the temperature and the boost factor we need to compare our data, not to 
the standard Maxwell distribution, but rather a Lorentz boosted Maxwell 
distribution. 

In principle this is a straight forward exercise, but it becomes complicated 
because the different components of the velocity couple through the Lorentz 
factor. Then, the Maxwell distribution of a Lorentz boosted thermal pop- 
ulation is not merely the Lorentz boost of the Maxwell distribution of the 
population at rest. 

Below in Eq. (A. 11) and Eq. (A. 12) I present the Maxwell distribution 
functions as function of the boost factor F, boost velocity u and temperature 



The standard Maxwell distribution for a population at rest in its most basic 
form can be written 



where dN is the number of particles per dvxdvydvz and N{T) is an overall 
normalisation factor. Going to spherical coordinates and integrating out the 
angle dependence it changes to 



T. 



A.l The standard relativistic distribution 




(A.l) 




(A.2) 



A.2 Boosting the Maxwell distribution 
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while the most convenient system for boosting the distribution is cyhndrical 
coordinates, where it can be written 

dN = 27rA^(T) exp ^-^^^^ v±dv±dvz . (A.3) 

When considering PDFs, from a numerical point of view, the most nat- 
ural variable to work in is not the normal velocity. The three velocity is 
bounded by the speed of light and the PDFs are squeezed towards c at high 
temperatures. Instead normally the four velocity vy is used, which is linear 
all the way from non relativistic to ultra relativistic velocities. The Maxwell 
distribution in terms of and 7 is 



dN = 47rA^(r)^^^2^ exp (^-^^ ) (A.4) 



and 



- 47rAr(T)^^-i|^^ exp (-1^) ^(.7) ■ (A.5) 

A.2 Boosting the Maxwell distribution 

To generalise the above distributions to those seen by observers moving with 
four velocity uT along the 2;-axis we need to Lorentz transform the vari- 
ables. The Lorentz transformation together with the inverse transformation 
between the two rest frames are 

J = Tj^l — uvz) 7 = r7'(l + Mf^) (A. 6) 

' ^ vz-^ (A.7) 



1 — UVz 1 -|- UV'^ 

^±- p.-,l^ ^ v'± = '"^ . , (A.8) 

r(i + uv'j r(i - uvz) 

where v± is a velocity component perpendicular to the boost direction and 
prime denots the boosted reference frame. To derive the Maxwell distribu- 
tion, as seen by an observer moving in the z-direction, we have to transform 
either Eq. (A.2) or Eq. (A.3) and reexpress it in terms of the new coordi- 
nates. The Maxwell distribution in cylindrical coordinates is best suited, 
since from Eq. (A. 6) we see that the transformation of 7 will pick up a de- 
pendence on v'^. Using Eqs. (A.7) and (A.8) to evaluate the Jacobian of the 



A.3 The boosted Maxwell velocity distribution 
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differentials and substituting the new variables into Eq. (A.3), the boosted 
Maxwell distribution in cylindrical coordinates may be found to be 

dN = 27rN(T) exp /^_ ry(l + - 1 \ < Mdv', . (A.9) 

In this form the distribution function cannot be compared directly with PDFs 
obtained from numerical data, since it is still two dimensional. We need to 
marginalise one of the two dimensions, to reduce it to a one dimensional 
PDF. 



A.3 The boosted Maxwell velocity distribution 

To find the velocity distribution we shift to spherical coordinates, setting 

v', = v' cos(e) v'^ = v' sm(e) , (A.IO) 

where 9 G [0,7r],w' G [0,1]. Inserting the new coordinates in Eq. (A.9) and 
integrating over angles after some algebra the boosted Maxwell distribution 
binned linearly in the velocity is 

,^ = 2.7V(T)Tl^y___ (A.11) 

where the temperature dependent integral has the limits a± = ^^^^^ , 

As mentioned above, when analysing particle data it is important to 
compute the PDFs with a linear behaviour from sub- to ultra relativistic 
velocities. Changing from dv' to d^v' we find the final result 

The integral in Eq. (A.12) may be simplified by repeated partial integration 

r e-Pd(5 _ -(i + r/3)2 + r(i + r/3) -2r2 i r 

JJTTtW' 673(1 + Tpf ~&r^J (l + T/3) 

(A.13) 

and everything reduces to an exponential integral, that depends on T. 

When analysing data I use IDL. It already contains a function to evaluate 
the exponential integral, and it is rather trivial to implement Eq. (A.12) into 
a computer program that, given a set of particles, evaluates the PDF, fits a 
boosted Maxwell distribution and finds the corresponding temperature and 
velocity. 



B. TRANSFORMATION OF TENSORS BETWEEN 
DIFFERENT METRICS 



In Chapter 2 it was argued that calculating variables in a local frame, re- 
taining at the same time global coordinates is the best approach for our nu- 
merical method. Methods used for special relativity may then be employed 
with minimal changes in arbitrary space times. 

In this appendix, I give the detailed transformation rules for vectors and 
two tensors. I consider the transformation between three different frames. 
The global star fixed coordinate system (SFCS) has the metric 

ds"^ = -a^dt^ + 7y {dx' + (3'dt) {dx^ + p^dt) . (B.l) 
The local laboratory (LOLA) frame has the metric 

ds'^ = -dt^ + -fijdx'dx^ , (B.2) 
while the pseudo fiducial observer (PFIDO) frame has the metric 

ds^ = -dp + J2 -^dx'dx^ . (B.3) 

In the case that the metric contains no off diagonal spatial components, the 
PFIDO frame is, in fact, the frame of a fiducial observer. In the worst case 
the PFIDO metric contains three non-trivial components. 
The three metrics are related by the relations 

{dx' + (3'dt) = dx' adt = di (B.4) 

^/^idx^ ^ dx^ di^di. (B.5) 

The differentials transform as contravariant vectors. The transformation laws 
for contravariant vectors may be found by multiplying with metrics and doing 
a bit of linear algebra. Tensors by definition transform as the product of the 
corresponding vectors and it is a straight forward, though tedious, exercise 
to find all the combinations. I have written them down here, since they 
are essential for the implementation of any physics; deriving how the local 
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variables are related to the global one. The following relations have been of 
interest when transforming to and from different frames: 
SFCS ^ LOLA frame: 
(vectors) 

(B.6) 
(B.7) 

(B.8) 

(B.9) 
(B.IO) 

SFCS ^ LOLA frame: 

(contravariant two-tensors) 
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(mixed type two-tensors) 
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a 

rpt ^ rpt 

'~ a ' 

Ti = a (fi - (^fn + (3^ ( r; - ^r] 



a J \ a 



/3« . 
3 ^ a 



(covariant two-tensors) 



Ttt 




Tti 




T 


— f 



(3^^ \ , , (3- 



— -^tj -r u:/-- -tit -I J-ij 

a \ a 



a 



B.ll) 
B.12) 

B.13) 



B.14) 
B.15) 
B.16) 

B.17) 

B.18) 

B.19) 
B.20) 
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LOLA ^ SFCS frame: 

(contravariant two-tensors) 

ftt ^ ^2ri.tt ^g_2l) 

T*^ = a (T*^ + (B.22) 

(mixed type two-tensors) 

t/ = T/ - p'T^ (B.24) 

T;^ = aT^ (B.25) 

= 1 [t; + /3^7;* - (t; + /3^T/)] (B.26) 

f; = Tj + (B.27) 

(covariant two-tensors) 

fu = ^ {Ttt - P'Ttj - P'Tu + P'P'Tij) (B.28) 

fti = - {Tti - 13%,) (B.29) 
a 

fij = Tij (B.30) 

LOLA frame ^ PFIDO frame: 

(vectors) 

[>* = [/* C/^ = ^[/^ (B.31) 

t/* = t/* C/^ = y/^tj' (B.32) 
(contravariant two-tensors) 

rj-itt rj^ti rpti 33) 



/7i> 

T^^' = ^ (B.34) 

(mixed type two-tensors) 

fl = fl fl = (B.35) 
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(covariant two-tensors) 



Ttt — Tft 

Tij = \/%i^yjjTij 



Tti — "V^ Ifii-T'ti 



(B.37) 
(B.38) 
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